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Abstract 

Since  2004,  AFIT  has  been  developing  a  trend- analysis  tool  to  assess  large  com¬ 
mercial  aircraft  infrared  (LCAIR)  signatures.  In  many  cases,  this  code  predicted 
signatures  to  within  10%  of  measured  data.  However,  other  results  indicated  that  the 
single-bounce,  specular-reflection  algorithm  being  used  failed  to  adequately  simulate 
interactions  between  aircraft  parts  where  either  the  specular  component  is  dominated 
by  diffuse  reflection  or  part-to-part  multiple-bounce  reflections  contribute  significantly 
to  the  signature;  discrepancies  greater  than  100%  were  observed.  This  research  in¬ 
corporates  Bi-Directional  Reflectance  Distribution  Functions  (BRDF’s)  and  multiple- 
bounce  calculations  into  the  LCAIR  model. 

A  physical  aircraft  model  was  constructed  from  aluminum,  and  measurements 
were  taken  before  and  after  a  surface  treatment  in  gloss  black  paint.  The  Sandford- 
Robertson  model  is  used  to  parameterize  the  BRDF’s  of  both  the  bare  aluminum  and 
gloss  black  paint.  Since  the  most  efficient  method  of  integrating  a  BRDF  depends 
upon  the  reflectance  distribution  of  the  aircraft  material,  the  sampling  resolution  of 
the  BRDF  integral  is  crucial  to  an  accurate  simulation.  Additionally,  care  is  taken 
to  ensure  that  the  integration  of  the  hemispherical  irradiance  onto  each  facet  of  the 
computational  model  is  sampled  at  a  sufficient  resolution  to  achieve  convergence  in 
the  solution. 

Simulations  in  the  mid-wave  infrared  (MWIR)  and  long-wave  infrared  (LWIR) 
bands  validate  both  the  previous  specular  reflectance  simplification  for  the  gloss  black 
simulations  and  the  failure  of  the  previous  algorithm  for  the  highly  reflective  bare 
aluminum.  The  necessity  of  considering  multiple  bounces  in  the  simulation  is  also 
demonstrated  amongst  part-to-part  reflections  near  the  wing  root,  where  three  or 
four  bounces  are  required  for  the  solution  to  converge. 


IV 


Finally,  three  scenarios  simulating  a  man-portable  air  defence  system  (MAN- 
PADS)  system  engaging  an  Airbus  A340-300  aircraft  landing  at  a  generic  airport  are 
performed.  The  infrared  signature  of  the  aircraft  is  contrasted  against  the  background 
over  the  detector’s  field  of  view.  An  irradiance  at  the  detector’s  optic  is  calculated 
for  the  MWIR  and  LWIR  bands.  A  maximum  noise  equivalent  irradiance  (NEI)  of 
5  x  10-5  W/m 2  is  found  to  be  required  to  ensure  that  the  MANPADS  is  able  to  detect 
and  track  the  aircraft.  An  NEI  larger  than  this  will  prevent  the  system  from  detecting 
the  aircraft. 
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A  Study  of 


Bi-Directional  Reflectance  Distribution  Functions 

and  Their  Effect  on 
Infrared  Signature  Models 

I.  Introduction 

1 . 1  Infrared  Spectrum 

The  infrared  (IR)  spectrum  is  a  subset  of  the  electromagnetic  spectrum  between 
visible  radiation  (light)  and  radio  waves  (millimeter  and  longer  wavelengths) 
(Figure  1.1).  It  has  been  used  extensively  in  the  scientific  community,  with  examples 
of  communications,  spectroscopy,  and  the  humble  remote  control. 

The  most  common  and  widespread  application  is  thermal  emission  and  imaging. 
Industries  including  home  security,  meteorology,  night  vision  and  thermal  heating  rely 
on  the  infrared  spectrum’s  physical  properties  for  their  operation.  Another  application 
of  infrared  energy  is  infrared  guided  (heat  seeking)  missiles.  These  missiles  exploit 
the  property  that  all  matter  above  0  Kelvin  emits  radiation.  This  energy  is  known  as 
Planckian  radiation,  and  for  most  objects  at  realizable  temperatures,  the  emission  is 
in  the  IR  spectrum. 
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Figure  1.1:  The  infrared  spectrum  is  a  subset  of  the  electro¬ 
magnetic  spectrum.  Wavelength  is  in  units  of  microns  (/im). 
Adapted  from  [4,  Figure  1.1]. 
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(a)  (b) 

Figure  1.2:  The  DHL  A330  aircraft  was  hit  by  an  SA-14  MAN- 
PADS  in  the  left  wing  and  lost  all  hydraulics.  Reproduced  from 
airdisaster.com. 

1.2  The  MANPADS  Threat 

Since  the  proliferation  of  thousands  of  man  portable  air  defence  systems  (MAN- 
PADS)  after  the  fall  of  the  Soviet  Union,  there  has  been  a  credible  threat  to  commer¬ 
cial  and  military  aviation  from  terrorist  and  issue  motivated  groups  (IMGs).  Several 
aircraft  have  been  shot  down  in  past  years,  and  whilst  these  events  have  occurred 
mainly  in  European  and  African  countries,  the  current  terrorist  threat  level  to  the 
United  States  and  other  allied  countries  is  high.  A  recent  incident  occurred  when  a 
DHL  Airbus  A330  cargo  aircraft  was  hit  by  an  SA-14  MANPADS  after  takeoff  from 
Baghdad  International  Airport  on  22  November  2003.  The  damage  was  so  severe  that 
the  aircraft  lost  all  three  hydraulic  systems,  and  the  aircrew  were  forced  to  control 
the  aircraft  with  only  differential  thrust.  This  was  the  only  time  an  aircraft  has  safely 
landed  with  all  hydraulic  systems  inoperable.  In  all  previous  attempts,  the  aircraft 
crashed  catastrophically  with  major  loss  of  life.  Figure  1.2  shows  the  damage  that 
the  DHL  aircraft  sustained  from  the  missile  hit. 

The  issue  is  recognized  by  governments  around  the  world  and  the  US  Depart¬ 
ment  of  Homeland  Security  has  expressed  a  desire  to  mandate  IR  countermeasures 
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be  installed  on  large  commercial  aircraft  in  the  US.  Legislation  has  been  introduced 
into  the  US  Congress  for  consideration,  (HR.  580,  S.311). 

1.3  Problem  Description 

This  research  topic  began  in  2004  with  Capt  Ruben  Martinez  who  created  a 
trend  analysis  tool,  named  LCAIR  (Large  Commercial  Aircraft  Infrared  Trend  Anal¬ 
ysis  Tool)  [17].  The  tool  was  developed  with  the  intent  to  quickly  calculate  IR  signa¬ 
ture  trends.  The  trends  were  hoped  to  enable  low  observability  properties  to  be  better 
designed  into  an  aircraft  during  the  design  stage,  rather  than  as  an  afterthought. 

Capt  Jonathan  Bortle,  in  2006,  validated  the  LCAIR  code  using  a  physical 
model,  where  calibrated  measurements  were  made  for  comparison  with  the  LCAIR 
prediction  for  the  same  parameters.  The  code  performed  well  and  was  over  90% 
accurate,  however,  it  failed  in  certain  situations  when  a  highly  reflective  surface  was 
used.  Additionally,  LCAIR  only  uses  a  single  bounce  specular  reflectance  algorithm, 
which  provides  some  other  limitations  that  will  be  discussed  in  Chapter  III. 

1.3.1  Purpose.  This  thesis  investigates  the  effect  of  Bi-Directional  Re¬ 
flectance  Distribution  Functions  (BRDFs)  when  incorporated  into  IR  signature  sim¬ 
ulations,  with  the  intent  to  analyze  and  improve  on  the  limitations  of  LCAIR  that 
Bortle  described  in  [3] .  Additionally,  several  algorithms  will  be  investigated  to  deter¬ 
mine  the  optimum  balance  between  simulation  accuracy  and  calculation  time. 

1.4.  Thesis  Overview 

The  thesis  is  divided  into  five  chapters.  Chapter  I  introduced  the  research  topic, 
the  reasons  for  it,  previous  research  and  goals  for  this  thesis.  Chapter  II  will  introduce 
or  refresh  the  reader  with  concepts  and  background  theory  of  physical  phenomena, 
mathematical  methods  and  models  that  are  applied  in  the  thesis. 
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Chapter  III  develops  the  method  and  algorithms  that  will  be  used  to  evaluate 
the  different  techniques  utilized  to  calculate  the  IR  signature.  A  detailed  discussion 
and  development  of  the  simulation  parameters  and  assumptions  is  also  included. 

Chapter  IV  provides  the  detailed  results  of  each  of  the  simulations.  Three 
different  simulations  are  performed; 

•  Initially,  analytical  simulations  are  performed  to  quantitatively  validate  the  ac¬ 
curacy  of  the  algorithms.  These  are  followed  by  a  simulation  which  reproduces 
the  properties  of  an  artificial  cavity  blackbody  source, 

•  The  second  simulation  investigates  the  limitations  of  the  LCAIR  algorithm,  by 
including  (BRDF)  reflections  and  multiple  bounces  into  the  calculations.  The 
effect  of  different  combinations  of  these  is  also  shown. 

•  The  simulations  conclude  with  a  scenario  of  a  MANPADS  engaging  an  Airbus 
A340-300  aircraft  whilst  landing  at  a  generic  airport.  Three  positions  are  chosen 
and  the  IR  signature  at  the  three  observation  angles  is  calculated  to  determine 
the  required  parameters  of  the  MANPADS  to  detect  and  lock-on  to  the  aircraft. 

Chapter  V  provides  a  summary  of  the  results,  recommendations  for  further  work 
and  lessons  learned  over  the  course  of  this  thesis. 

Several  appendices  are  also  included  which  provide  information  that  helps  the 
reader  with  mathematical  methods  and  other  techniques  that  are  used  in  the  simula¬ 
tions. 

1.5  Literature  Review 

A  comprehensive  literature  review  was  performed  whilst  searching  for  IR  sim¬ 
ulation  techniques  and  propagation  algorithms.  Whilst  many  papers  are  available  in 
the  open  literature,  their  content  is  concerned  with  either  system  level  design  [7,8], 
or  focusses  on  specific  properties  of  BRDFs  [2,24], 
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Specific  code  and  algorithms  of  these  IR  simulation  tools  were  not  found,  as 
they  are  either  classified  or  commercial-in-confidence.  For  this  reason,  the  computer 
graphics  algorithms  found  in  [5]  and  [10]  were  used  to  improve  the  LCAIR  propagation 
algorithm  of  [3]  and  [17]. 

Several  papers  were  found  late  in  the  thesis  effort,  [14,16,21,22],  These  pa¬ 
pers  provide  a  good  introduction  into  IR  signature  modelling  and  analyze  important 
aspects  of  simulation  modelling,  trend  analysis  and  system  parameters.  Another  pa¬ 
per,  [13],  discusses  signature  analysis  in  the  context  of  aircraft  conceptual  design.  The 
content  of  these  papers  was  not  included  in  this  thesis  due  to  limited  time,  however, 
they  would  provide  a  good  starting  point  for  the  next  IR  signature  student. 


5 


II.  Background 


This  chapter  provides  applicable  background  theory  and  information  for  the 
reader  that  is  useful  in  understanding  the  concepts  and  methods  described  in 
Chapter  III. 

2. 1  Blackbody  Radiation 

It  is  assumed  that  the  reader  is  comfortable  with  the  concepts  of  blackbody 
radiation,  emissivity  and  reflectance.  The  relevant  equations  will  be  included  in  this 
chapter  for  reference  in  subsequent  chapters.  However,  if  the  reader  requires  extra 
information  relating  to  the  derivation  of  the  equations,  and  concepts,  they  are  directed 


to  [3,4,17]. 


A  blackbody  is  a  perfect  emitter  of  Planckian  radiation.  Kirchhoff ’s  law  states 
that  for  a  body  to  be  in  thermal  equilibrium  the  absorbed  equals  the  emitted  energy. 
Thus  the  incident  flux,  coincident,  equals  the  sum  of  the  absorbed,  reflected,  and  trans¬ 
mitted  flux.  If  the  body  is  opaque,  the  transmitted  flux  equals  zero,  and  then  the 
re-emitted  and  reflected  flux  equals  the  incident  flux: 


(2.1) 


A  perfect  blackbody  is  one  that  absorbs  all  incident  radiation.  To  remain  in  ther¬ 
mal  equilibrium,  the  body  must  re-emit  that  radiation.  The  spectrum  of  the  emitted 
radiation  is  dependent  upon  the  temperature  of  the  material.  Equation  (2.2)  shows 
the  expression  for  spectral  radiance,  Le(A),  of  a  perfect  blackbody  of  temperature, 
T,  in  Kelvin,  where  h  is  Plancks’s  constant,  c  is  the  speed  of  light,  k  is  Boltzmann’s 
constant,  and  A  is  the  wavelength  of  interest.  Figure  2.1  shows  the  emitted  spectral 
radiance  of  a  blackbody  (BB)  plotted  at  several  temperatures. 


Le  BB  Self-Emitted 
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Figure  2.1:  Blackbody  radiance  is  a  function  of  both  temper¬ 
ature  and  wavelength.  These  theoretical  curves  will  be  scaled 
by  emissivity,  s{6),  to  produce  greybody  emission. 

2.1.1  Reflectance.  Reflectance,  p{6 ),  is  defined  as  the  ratio  between  incident 
flux  and  reflected  flux  at  a  given  incident  elevation  angle,  0,  from  the  surface  normal. 
The  remaining  flux  is  absorbed. 


P(0) 


^Reflected 

^Incident 


(2.3) 


2.1.2  Emissivity.  Recall  from  Equation  (2.1)  that  what  is  not  reflected  is 
absorbed.  The  absorptance  factor  is  the  percentage  of  flux  that  is  absorbed  and  re¬ 
emitted,  and  for  an  opaque  surface  in  thermal  equilibrium,  is  defined  as  the  emissivity, 


m  =  i  -  m- 


(2.4) 


2.2  Radiometry 

Radiometry  is  the  mathematics  of  calculating  energy  (electromagnetic  flux) 
transfer  from  a  source  to  a  detector.  Exitant  flux  from  a  source,  incident  flux  onto 
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Table  2.1:  Integral  Forms  of  Radiance.  The  e  sub¬ 
script  implies  power  units  of  Watts,  not  photon  units. 
As  is  the  area  of  the  source  and  is  the  solid  angle 
subtended  by  the  detector. 


Quantity 

Variable 

Expression 

Units 

Intensity 

Ie 

f  Le  cos (9s)dAs 

As 

Exitance 

Me 

J*  cos {0 

^ d 

Flux 

$e 

J*  J  cos {0 

^ d  As 

m 

Irradiance 

Ee 

d^e 

dAd 

IS 

a  surface  and  collected  flux  by  a  detector  can  easily  be  calculated  and  converted  us¬ 
ing  radiometric  quantities.  The  most  basic  quantity  is  radiance,  and  carries  units  of 
Watts  per  steradian  per  unit  area.  Other  quantities  can  be  gained  from  radiance  and 
are  expressed  in  Table  2.1;  where  6S  is  the  angle  from  the  observation  angle  to  the 
surface  normal,  dAs  is  the  differential  area  of  the  source,  dQ,j  is  the  differential  solid 
angle  of  the  detector,  and  <Fe  is  the  integrated  flux  through  the  detector. 

2.3  Bi-Directional  Reflectance  Distribution  Function 

The  reflectance  property,  p(0),  represents  the  proportion  of  energy  that  is  re¬ 
flected  off  the  surface  for  a  given  incident  angle,  6,  however,  there  is  no  other  infor¬ 
mation  available  with  respect  to  the  reflection  angle.  The  specular  or  diffuse  nature 
of  the  surface  is  unknown.  In  practice,  energy  is  not  perfectly  specularly  reflected 
nor  scattered  equally  in  all  directions.  The  BRDF  provides  a  method  of  describing 
reflectance  as  a  function  of  incident  (source)  and  reflected  (observed)  angles  and  wave¬ 
length.  The  BRDF  is  defined  as  the  differential  ratio  of  the  reflected  radiance,  Lr,  to 
the  incident  irradiance,  £).  Nicodemus  et  al.  introduced  the  concept  of  the  BRDF 
in  [20]  which  contains  a  more  rigorous  treatment  of  the  definition,  assumptions  and 
applications.  Although  differential  quantities  cannot  be  directly  measured,  as  a  finite 
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Figure  2.2:  Hemispherical  geometry  of  the  BRDF,  defining 

the  elevation  and  azimuth  angles  for  both  the  source  and  obser¬ 
vation  vectors. 

solid  angle  of  energy  must  be  collected,  they  can  be  calculated  from  the  empirical 
data  of  each  surface.  As  defined  above,  the  BRDF,  with  units  of  inverse  steradians, 
is  shown  in  Equation  (2.5)  with  angles  as  defined  in  Figure  2.2; 


dLr(6r ,  (f)r  ) 


dEi(9i,<j>i ) 


w 

/ 

~W~ 

'  1  ' 

sr  —  m2 

m 2 

sr 

(2.5) 


where  9t  is  the  incident  elevation  angle  measured  from  the  surface  normal,  ch  is  the  in¬ 
cident  azimuthal  angle  measured  from  the  x-axis  or  some  reference,  6r  is  the  reflected 
elevation  angle  again  measured  from  the  surface  normal,  and  (j)r  is  the  reflected  az¬ 
imuthal  angle  measured  from  the  same  reference  as  </>*. 
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2.3.1  Dimensionality  Reduction.  The  full  BRDF  is  a  hi-fidelity  representa¬ 
tion  of  a  surface’s  reflectance  properties;  however,  it  is  very  complex  and  requires  a 
large  amount  of  processing  power  and  memory,  which  is  not  needed  for  some  appli¬ 
cations.  The  domain  of  the  BRDF  can  be  significantly  reduced  using  the  theory  of 
reciprocity,  an  assumption  of  surface  isotropy  and  wavelength  simplifications. 

2. 3. 1.1  Reciprocity.  The  theory  of  reciprocity  states  that  the  path 
of  light  is  independent  of  the  direction  of  propagation.  Reversing  the  direction  of 
propagation  does  not  alter  the  reflectance  properties.  The  BRDF  domain  is  reduced 
by  using  the  symmetry  of  the  variables, 


frifiii  Pit  Pri  fri^ri  Pri  Pii  ^)-  (2.6) 

2. 3. 1.2  Wavelength.  If  the  BRDF  can  be  assumed  to  not  change 
significantly  over  a  wavelength  band,  the  BRDF  can  be  set  constant  across  the  band. 
This  assumption  breaks  down  as  the  spectral  width  of  the  band  increases  where  the 
majority  of  materials  become  more  specular  as  the  wavelength  increases.  This  effect 
will  be  clearly  seen  later,  in  Figure  2.7,  where  the  diffuse  component  decreases  and 
specular  component  increases  as  the  wavelength  increases. 

2.3. 1.3  Surface  Isotropy.  Surfaces  can  be  considered  isotropic  if  their 
reflective  properties  are  independent  of  the  azimuthal  incidence  angle,  pi]  most  com¬ 
mon  surface  types  fit  this  definition.  Two  exceptions  to  this  assumption  are  polar¬ 
ized  materials  and  materials  with  surface  striae  (e.g.,  gratings).  These  materials  will 
naturally  exhibit  different  reflectance  characteristics  as  a  function  of  pi.  Under  the 
isotropic  assumption,  the  strict  dependence  upon  pi  and  pr  with  reference  to  the 
x-axis,  as  shown  in  Figure  2.2,  is  removed  and  only  the  delta  (A p  =  pi  —  pr)  is 
required.  Applying  these  assumptions  to  Equation  (2.5)  gives  the  simplified  BRDF 
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with  three  dimensions  ($,,  0r  and  A </>), 


fr(0i,  0r,  A 0)  =  /r(0i,  &,  0r,  0r,  A).  (2.7) 

Integral  Forms  of  the  BRDF.  Equation  (2.5)  can  be  converted  to 
many  quantities  by  integrating  across  one  or  more  domains  of  the  BRDF.  The  integral 
forms  provide  information  about  the  material  that  may  not  be  intuitive  by  piecewise 
evaluation  of  the  BRDF.  A  few  common  quantities  that  will  be  used  later  are  defined 
below  and  a  full  list  of  possible  integrated  quantities  can  be  found  in  [20,  Tables  2  &  3]. 

2.3.2. 1  Directional  Hemispherical  Reflectance.  The  Directional  Hemi¬ 
spherical  Reflectance  (DHR)  describes  the  percentage  of  the  energy  reflected  into  the 
entire  hemisphere  from  a  given  incident  angle  [20,  Table  3,  Equation  (3.1)].  This 
function  is  derived  from  the  integral  of  the  BRDF  over  all  exitant  angles  as  shown 
in  Equation  (2.8),  where  the  </>*  dependence  is  dropped  under  the  previous  surface 
isotropy  assumption.  However,  it  is  important  to  keep  the  </>;  reference  as  a  baseline 
for  the  propagation  algorithms  described  in  section  3.5.  The  DHR  must  also  equal 
the  reflectance,  p(0),  as  reflectance  is  the  ratio  of  reflected  flux  to  incident  flux.  In¬ 
tegrating  the  BRDF  over  all  reflected  angles  accounts  for  all  reflected  energy  (flux), 
and  is  thus  a  ratio  total  energy. 

27 r  tt/2 

DHR{0i)  =  piOf)  =  J  j  fr(0i,  fa,  0r,  4>r)  cos(6*r)  sm(9r)d9rd(j)r  (2.8) 
o  o 

2. 3. 2. 2  Hemispherical  Directional  Reflectance.  Conversely,  the  Hemi¬ 
spherical  Directional  Reflectance  (HDR)  describes  the  energy  that  is  reflected  from 
the  whole  hemisphere  into  a  specific  direction  [20,  Table  3,  Equation  3.6].  The  HDR 
is  obtained  by  integrating  the  BRDF  over  the  incident  angle  domains  as  expressed  in 
Equation  (2.9)  with  the  azimuthal  dependence  removed  as  developed  before. 
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(a) 


(b) 


Figure  2.3:  (a)  A  purely  specular  reflection  does  change  the 

distribution  of  energy  from  incident  to  reflected  directions. 

(b)  Semi-specular  reflections  are  almost  specular  where  the  in¬ 
cident  distribution  spreads  slightly  and  can  be  described  with 
some  lobe  width. 


HDR(6r) 


27 r  tt/2 


fr(6i,  (pi,  0r,  4>r)  cos {Of)  sin (d^ddidfa 


0  0 


(2.9) 


2.3.3  Surface  Reflectance  Distributions.  Using  the  BRDF,  the  reflectance 
of  a  surface  can  now  be  described  in  detail  as  a  function  of  incident  and  reflected 
angles  to  represent  the  specific  reflectance  distribution  for  that  material  type.  Two 
trivial  cases  are  diffuse  and  specular  surfaces.  Other  cases  are  represented  by  the 
distribution  of  the  reflectance  properties  which  are  discussed  below. 


2.3.3. 1  Specular  Reflectance.  Specular  reflection  is  perhaps  the  most 
common  reflectance  type.  It  is  the  major  assumption  in  ray  optics  and  is  the  basis 
of  the  well  known  Law  of  Reflection,  6i  =  6r.  Surfaces  which  are  optically  smooth  (a 
material  with  surface  roughness  appreciably  smaller  than  the  wavelength  of  interest) 
can  be  considered  specular.  The  Law  of  Reflection  requires  that  the  angle  of  incidence 
equals  the  angle  of  reflection;  this  is  the  two  dimensional  (2D)  definition.  In  three 
dimensions  (3D),  the  elevation  angles  for  the  incident  and  reflected  rays  must  equal 
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Figure  2.4:  The  reflected  radiance  distribution  is  constant, 

fr  =  1  / 7r  steradians,  in  all  directions  for  the  perfectly  diffuse 
reflectance  distribution. 


and  the  azimuthal  angles  must  differ  by  7 r  radians.  Applying  these  conditions  to  the 


BRDF,  yields  Equation  (2.10),  and  is  depicted  in  Figure  2.3  (a).  The  specular  case 
is  special  and  usually  requires  a  separate  algorithm  to  evaluate  the  equality. 


(2.10) 


A  relaxation  to  strict  equality  conditions  of  the  specular  case  is  known  as  the  semi- 
specular  case.  Semi-specular  is  when  the  surface  is  almost  specular,  however,  there 
is  some  width  to  the  reflected  solid  angle,  producing  a  specular  lobe  instead  of  a 
perfectly  reflected,  collimated  beam.  This  effect  is  chiefly  caused  by  the  surface  not 
being  completely  optically  smooth  and  is  depicted  in  Figure  2.3  (b). 


2. 3. 3. 2  Diffuse  Reflectance.  The  opposite  of  the  specular  case  is  where 


the  energy  is  scattered  equally  in  every  direction  in  the  hemisphere,  as  depicted  in 
Figure  2.4,  and  is  defined  as  diffuse  reflectance.  Some  common  diffuse  materials  are 
matte  paints,  carpet  and  video  projection  screens.  As  the  exitant  radiance  from  the 
surface  is  constant  in  every  direction,  so  too  is  the  BRDF.  Recall  that  the  DHR  is 
the  integral  of  the  BRDF  over  all  exitant  angles  and  is  equal  to  the  reflectance  at 
that  incident  elevation  angle.  Therefore,  the  diffuse  BRDF,  when  integrated  over  the 
exitant  hemisphere,  must  equal  p(9i).  This  leads  to  the  value  of  1/7T  for  the  diffuse 
BRDF  and  is  scaled  by  p{0f), 
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Figure  2.5:  The  fractional  specularity  reflectance  model  is  a 

combination  of  a  diffuse  and  semi-specular  surface. 


fr(9i,0r,  A(f>) 

II 

(2.11) 

27T  7t/2 

r  r 

DHR(9i,  fa) 

=  /  /  fr(9i,  fi,  9r,  (fir)  cos(9r)  sin(9r)d9rdfr 

(2.12) 

0  0 

27 r  7r/ 2 

cos(6)r)  sin (Or)dOrd<frr  =  p{9i).  (2-13) 


o  o 

2. 3. 3. 3  Fractional  Specularity  Reflectance  Model  .  In  reality,  few  sur¬ 
faces  fit  the  definitions  of  perfectly  diffuse  or  specular,  thus,  simple  analytical  de¬ 
scriptions  of  the  BRDF  are  not  sufficient  to  represent  the  reflectance  distribution 
throughout  the  hemisphere.  For  surfaces  that  are  nearly  specular,  but  exhibit  some 
diffuse  reflectance  (or  vice-versa),  a  fractional  specularity  model  can  be  defined.  The 
model  can  be  defined  in  two  ways;  either  the  hemisphere  is  separated  into  diffuse 
and  specular  regions  in  a  piecewise  approach,  or  a  complete  hemispherical  diffuse 
component  is  summed  with  the  specular  component  to  represent  the  distribution;  an 
example  of  the  second  method  is  shown  in  Figure  2.5. 

2. 3. 3. 4  Parameterized  BRDF  Models.  Sections  2.3.3. 1  through  2. 3. 3. 3 
discussed  theoretical  models  to  represent  reflectance  distributions  of  surfaces.  How¬ 
ever,  real  surfaces  do  not  exhibit  uniform  distributions  like  the  diffuse  case  or  the  delta 
function  of  the  specular  case;  thus,  their  empirical  reflected  radiance  measurements 
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will  vary  with  elevation  and  azimuth  angles.  Materials  (especially  specular  ones)  can 
exhibit  large  gradients  in  the  BRDF  function,  which  requires  a  high  sampling  rate 
in  the  empirical  measurements.  This  ensures  the  reflectance  distribution  is  recorded 
properly,  however,  this  process  produces  large  amounts  of  data. 

For  practical  uses,  it  may  not  be  feasible  to  store  the  entire  array  of  measured 
data  for  recall.  Additionally,  a  BRDF  value  for  a  point  which  lies  between  sample 
measurements  may  be  required.  In  these  two  cases,  it  is  much  simpler  to  develop 
a  model  which  parameterizes  the  BRDF  surface  allowing  for  fast  evaluation  of  the 
required  value,  independent  of  initial  measurement  locations.  The  functions  to  be  fit 
to  the  empirical  data  depend  very  much  upon  the  empirical  BRDF  shape;  a  single 
model  will  not  give  good  results  for  all  surface  reflectance  types.  Obviously,  the  more 
empirical  measurements  and  the  more  degrees  of  freedom  in  the  functions,  the  more 
accurate  the  parametrization  will  become. 

There  are  many  BRDF  models  available  to  parameterize  empirical  BRDF  mea¬ 
surements  that  can  be  found  in  literature.  Sandford- Robertson  [23],  Maxwell-Beard  [18], 
and  He-Torrance-Sillion  [12],  are  just  three.  This  thesis  will  concentrate  on  the 
Sandford-Robertson  model,  as  it  is  predominantly  used  for  IR  BRDF’s  on  common 
aircraft  surfaces.  The  Air  Force  Research  Laboratory  (AFRL)  sponsored  the  devel¬ 
opment  of  the  Sandford-Robertson  model. 

2. 3-4  Sandford-Robertson  Reflectance  Model.  One  of  the  common  reflectance 
models  used  in  aircraft  IR  signature  simulations  is  the  Sandford-Robertson  Model. 
Brian  Sandford  and  David  Robertson  published  the  model  in  1985  [23].  The  model 
is  a  fractional  specularity  model  and  is  divided  into  two  components,  specular  and 
diffuse. 

The  model  is  described  by  four  parameters,  which  represent  the  physical  prop¬ 
erties  of  the  material; 

•  Diffuse  spectral  reflectance,  pp( A), 
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•  Spectral  emissivity,  e(A), 

•  Grazing  angle  reflectivity,  6;  and 

•  Width  of  specular  lobe,  e. 

The  functions  to  parameterize  the  empirical  data  are  separated  into  diffuse  and  spec¬ 
ular  components,  where  the  above  parameters  are  chosen  to  fit  both  a  diffuse  floor 
and  a  specular  lobe  to  the  measured  data.  The  total  BRDF,  fr,  is  the  sum  of  the 
diffuse,  fd,  and  specular,  fs,  components, 

fr  =  fd  +  fs-  (2-14) 


2.3.4- 1  Diffuse  Floor.  The  diffuse  floor  is  assumed  to  be  symmetric  in 
azimuth  as  a  result  from  subsurface  scattering  and  surface  roughness.  It  is  calculated 
by; 


,  =  i  g(Or)pD(X)g(0i) 
Jd  vr  [( G(b )]2 

where  the  g(9r),g(6i )  and  G(b)  are  defined  by, 


(2.15) 


g(9r) 

g{0i) 

Gib) 


l 

1  +  b2  tan  2{6r) 

1 


1  +  b2  tan 2(0i) 


1 

1  -  b2 


1 


b2 

1-b2 


logf) 


(2.16) 

(2.17) 

(2.18) 


For  a  more  detailed  explanation  of  the  components,  see  [23]. 


2. 3. 4-2  Specular  Lobe.  The  specular  reflectance  in  the  Sandford- 
Robertson  model  is  a  modified  version  of  the  Trowbridge  and  Reitz  model  [25].  The 
lobe  is  calculated  by  Equation  (2.19)  and  is  normalized  so  the  hemispherical  integral 
of  the  BRDF  is  equal  to  the  DHR  measurements.  Equations  (2.20)  through  (2.24) 
define  the  components  that  calculate  fs.  Again,  the  reader  is  referred  to  [23]  for  a 
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detailed  derivation  and  further  information  relating  to  the  Sandford  Robertson  BRDF 
Model. 


fa 


h(a)  1 
H(9i)  cos (9r) 


d(9r ) 


(2.19) 


Ps{\  @i) 


e(\,9i) 


h(a ) 


H{9i) 


d(9r) 


1  -  pd{\9i)  -  e{\ ,9i) 


£(A) 


Gib) 


1 


[e2  cos2  («)  +  sin2  (a)] ' 


1 

2^ 


(1  —  e2)  cos (9i)  + 


[2e2  +  (1  —  e2)2  cos2(6)i)] 
y/{l  —  e2)2  cos2(6,j)  +  4e2 


1  +  b2  tan 2(^r) 


(2.20) 

(2.21) 

(2.22) 

(2.23) 

(2.24) 


2. 3. 4-3  Glint  Vector.  The  majority  of  figures  illustrating  the  model’s 
goodness  of  fit  to  empirical  data  only  show  the  in-plane  reflectance  (</y  =  ±  7r). 

Figure  2.9  is  an  example  of  the  representation  style.  However,  values  are  required 
for  the  out-of-plane  angles  during  the  simulation.  The  glint  vector,  g ,  represents  the 
distance  from  the  specular  reflection  vector  and  is  specifically  defined  as  the  bisector  of 
the  incident  vector  and  the  reflected  vector.  It  is  shown  in  Figure  2.6  and  is  calculated 
by  Equation  (2.25).  The  glint  angle,  a,  is  then  calculated  by  Equation  (2.26)  for  use 
in  Equation  (2.22). 


9 

a 


(»+») 

V^2(l  +  6  •  s) 
cos-1(<7  •  h) 


(2.25) 

(2.26) 


The  Sandford-Robertson  model  parameterizes  with  the  assumption  that  the  surface 
is  predominantly  specular  with  a  diffuse  floor;  this  is  why  the  glint  vector  and  a 
angle  are  used  to  provide  a  larger  attenuation  in  BRDF  value  the  further  away  the 
observation  vector  is  from  the  specular  condition,  in  both  <j)r  and  9r. 
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Figure  2.6:  Sketch  of  the  glint  vector,  g ,  and  a  angle,  ref¬ 

erenced  to  the  source,  observation  and  normal  vectors.  The 
magnitude  of  the  a  angle  is  a  measure  of  the  angular  separation 
between  the  observation  and  specular  reflection  vectors.  Repro¬ 
duced  from  [23]. 

2.3.5  Bare  Aluminum  Material.  The  first  surface  treatment  of  interest 
is  bare  aluminum.  Bortle  measured  the  reflectance  data  using  DHR  measurements 
of  bare  aluminum.  The  results  are  included  in  Figure  2.7,  reproduced  from  [3]  for 
completeness.  These  spectral  curves  were  then  averaged  over  two  bands.  The  mid 
wave  infrared  (MWIR)  band  spans  3  —  5 fim,  and  the  long  wave  (LWIR)  band  covers 
8  —  12 fim.  Figure  2.8  provides  average  reflectance  data  as  a  function  of  incident 
elevation.  The  data  points  were  then  parameterized  with  polynomial  functions  to 
provide  the  p(6i)  functions  in  the  calculation  of  emissivity  and  self  emitted  radiance 
of  bare  aluminum  surfaces. 

BRDF  data  was  taken  with  AFRL’s  scatterometer  to  provide  an  empirical 
BRDF.  These  empirical  results  were  then  parameterized  by  the  Sandford- Robertson 
model;  both  sets  of  data  are  shown  in  Figure  2.9  with  parameters  as  listed  in  Ta¬ 
ble  2.2.  An  attribute  to  note  of  the  aluminum  surface  is  the  diffuse  component,  which 
is  considerable  at  3 pm,  but  decreases  to  almost  zero  at  12 pm.  Additionally,  as  the 
diffuse  component  is  reducing,  the  specular  component  increases  to  keep  the  total 
reflectance  relatively  constant  across  the  two  wave  bands. 
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Table  2.2:  Sandford-Robertson  BRDF  model  pa¬ 

rameters  for  the  bare  aluminum  surface,  reproduced 
from  [3].  A  detailed  approach  to  the  selection  and 
parameterizations  process  can  also  be  found  in  [3]. 


Spectral  Region 

e(A) 

Pd  (A) 

b 

e 

MWIR  (3  —  5 pm) 
LWIR  (8  -  12 pm) 

0.110 

0.080 

0.020 

0.020 

0.050 

0.050 

0.0042 

0.0042 

The  term  diffuse  component  is  misleading.  The  previous  definition  of  diffuse 
indicated  that  it  was  distributed  across  the  whole  hemisphere;  however,  the  DHR 
diffuse  measurement  strictly  means  non-specular.  When  the  DHR  measurements  are 
made,  a  small  baffle  is  inserted  into  the  reflectometer  to  block  the  specular  component 
to  measure  the  diffuse  component.  This  baffle  only  blocks  a  very  small  solid  angle, 
thus,  energy  that  is  still  in  the  specular  lobe  is  recorded  in  the  diffuse  section  of 
the  DHR  measurements.  Figure  2.9  shows  the  specularity  of  the  aluminum  material, 
it  just  happens  that  the  specular  lobe  is  wide  enough  to  extend  into  the  diffuse 
measurement  solid  angle. 

2.3.6  Gloss  Black  Painted  Aluminum  Material.  Similarly  to  the  bare  alu¬ 
minum  case,  DHR  (Figure  2.10),  p(6)  (Figure  2.11)  and  BRDF  (Figure  2.12)  mea¬ 
surements  were  taken  with  the  gloss  black  painted  aluminum  surface  and  the  BRDF 
was  parameterized  using  the  Sandford-Robertson  model  (Table  2.3). 

An  attribute  of  the  gloss  black  paint  to  note  is  the  magnitude  of  the  specular 
peak  of  the  BRDF ;  it  is  approximately  two  orders  of  magnitude  lower  than  that  of  the 
bare  aluminum  surface.  This  can  be  seen  in  the  DHR  data  also,  where  the  specular 
reflectance  of  the  gloss  paint  is  W1.3  compared  with  the  bare  aluminum  surface’s 
reflectance  of  ~0.8.  The  effect  of  this  will  be  shown  later  in  Chapter  IV. 
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Figure  2.7:  Directional  Hemispherical  Reflectance  measurements  for  the  bare  alu¬ 
minum  surface,  reproduced  from  [3,  Figure  3.8]. 


Table  2.3:  Sandford- Robertson  BRDF  model  pa¬ 

rameters  for  the  gloss  black  paint  surface,  reproduced 
from  [3].  A  detailed  approach  to  the  selection  and 
parameterizations  process  can  also  be  found  in  [3]. 


Spectral  Region 

e(A) 

Pd(  A) 

b 

e 

MWIR  (3  —  5/um) 
LWIR  (8  —  12/j.m) 

0.890 

0.880 

0.001 

0.00075 

0.160 

0.160 

0.005 

0.005 
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Figure  2.8:  p{6)  curves  for  bare  aluminum  in  both  MWIR  and  LWIR  bands.  The  y- 

axis  limits  are  set  to  [0  1]  as  a  reference  to  show  the  difference  between  wave  bands  and 
material  types.  For  images  showing  goodness  of  fit  for  the  reflectance  parametrization, 
see  [3,  Figure  3.9  (b)]. 
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Figure  2.9:  Empirical  and  parameterized  BRDF  for  the  bare  aluminum  surface, 

reproduced  from  [3,  Figure  3.22  (b)]. 
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Figure  2.10:  Directional  Hemispherical  Reflectance  measurements  for  the  gloss 

black  surface,  reproduced  from  [3,  Figure  3.7]. 
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Figure  2.11:  p{9)  curves  for  gloss  black  paint  in  both  MWIR  and  LWIR  bands. 

The  y-axis  limits  are  set  to  [0  1]  as  a  reference  to  show  the  difference  between  wave 
bands  and  material  types.  For  images  showing  goodness  of  fit  for  the  reflectance 
parameterizations,  see  [3,  Figure  3.9  (a)]. 
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Figure  2.12:  Empirical  and  parameterized  BRDF  for  the  gloss  black  surface,  repro¬ 
duced  from  [3,  Figure  3.22  (a)]. 
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Figure  2.13:  The  spectral  transmission  of  IR  energy  varies 

with  wavelength  as  different  molecules  absorb  the  energy  at  var¬ 
ious  positions  in  the  spectrum.  Reproduced  from  [19,  Figure 
7-1.10], 

2.4  Atmospheric  Transmission 

IR  energy  that  propagates  through  the  atmosphere  experiences  absorption  and 
scattering  due  to  the  molecular  and  aerosol  constituents.  Water  vapor,  carbon  diox¬ 
ide,  ozone  and  other  minority  gases  all  spectrally  absorb  radiation,  and  re-emit  the 
energy  as  blackbody  radiation  at  the  ambient  temperature.  Absorption  occurs  spec¬ 
trally  as  the  IR  wavelength  matches  the  various  electronic,  vibrational  or  rotational 
transitions  for  each  molecule.  Thus,  the  absorption  spectrum  is  very  dependent  upon 
wavelength,  [19,26],  and  is  shown  in  Figure  2.13. 

When  designing  a  detection  system,  it  is  important  to  understand  and  exploit 
these  spectral  regions  of  high  transmittance,  as  designing  a  detection  system  over  a 
band  where  there  is  high  absorption  will  significantly  reduce  the  system’s  detection 
range.  For  this  reason,  systems  are  designed  to  detect  in  the  short  wave  IR  (SWIR) 
(1  —  3 pm),  MWIR  (3  —  5 pm)  or  LWIR  (8  —  12 pm)  regions. 
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Figure  2.14:  An  example  transmission  spectrum  for  the  (a) 

MWIR  and  (b)  LWIR  bands  is  shown.  Note  the  transmission 
windows  where  IR  radiation  is  not  absorbed.  The  path  length  of 
this  simulation  was  set  to  infinity  at  45°  elevation  on  a  cloudless 
summer  day. 


Software  programs  such  as  MODTRAN  predict  the  spectral  transmission  as  a 
function  of  geographical  position,  time  of  day,  weather  conditions,  and  path  position 
and  length.  One  such  simulation  in  the  MWIR  and  LWIR  bands  is  included  in 
Figure  2.14  for  reference. 


2.5  Monte  Carlo  Integration 

Monte  Carlo  integration  is  a  technique  used  to  numerically  evaluate  integrals, 
and  is  more  efficient  than  standard  numerical  quadrature  techniques  for  multi- dimensional 
integrals. 


2.5.1  Numerical  Quadrature.  To  evaluate  the  integral,  /,  of  f(x)  over  the 
domain  x  €  [o,  b ]  analytically,  the  expression  is, 


I  = 


(2.27) 
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Figure  2.15:  Deterministic  numerical  quadrature  where  the 

function  is  sampled  at  evenly  spaced  intervals  and  approximated 
by  a  rectangle  across  the  interval. 

and  solve  by  integrating  the  function  f(x)  and  applying  the  limits.  However,  this 
approach  becomes  infeasible  for  multidimensional  functions  and  ones  which  cannot  be 
integrated  analytically.  Numerical  integration  is  performed  by  quadrature,  where  the 
continuous  function  is  approximated  by  rectangles  in  one  dimension,  solid  rectangles 
in  two  dimensions,  and  so  on. 

Each  domain  must  be  divided  into  N  usually  equal  sized  regions,  and  the  value 
of  f(x)  is  evaluated  at  each  x%  point.  Figure  2.15  illustrates  the  one  dimensional 
case  which  can  be  extrapolated  into  further  dimensions  as  required.  Each  value  /(x*) 
is  then  scaled  by  a  weighting  factor,  Wj,  and  summed  to  evaluate  the  integral.  For 
uniform  divisions,  the  weighting  factor  is  Wi  =  (b  —  a)/N  and  is  constant  V  i.  The 
estimate,  /,  of  Equation  (2.27)  then  becomes 

7  =  f»i  =  f  /(^~a)  (2.28) 

2=1  2=1 

which  calculates  the  mean  value  of  f(x)  over  the  domain  and  multiplies  by  the  width  of 
the  domain.  Intuitively,  as  N  increases,  I  becomes  a  closer  estimate  of  the  analytical 
result,  /.  The  distance  between  /  and  I  can  be  measured  by  the  variance  in  the 
solution,  a2. 
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To  extend  numerical  quadrature  into  d  dimensions,  Nd  samples  of  the  function 
are  required. 


2.5.2  Monte  Carlo  Method.  Monte  Carlo  integration  is  similar  to  the  numer¬ 
ical  case,  the  one  difference  is  that  the  samples,  ay,  are  taken  from  a  random  variable, 
not  uniform  divisions.  The  advantage  of  this  method  is  its  efficiency  in  evaluating 
multi-dimensional  integrals  where  some  number  of  samples  (n  <  Nd)  can  achieve  a 
result  where  I  is  a  suitable  estimate  of  I. 


The  estimate  for  Monte  Carlo  integration  is  very  similar  to  Equation  (2.28), 
except  the  weighting  factor  differs.  The  (6  —  a)  term  is  replaced  with  l/p(ay)  in 
Equation  (2.29),  where  p(ay)  is  the  probability  density  function  (PDF)  of  ay  occurring 
in  the  integral  domain. 


1  N 

;4e 


f(Xi) 


(2.29) 


N  “1 '  P{Xi) 

This  substitution  remains  consistent  with  Equation  (2.27)  as  the  integral  of  p(ay)  over 
the  uniformly  sampled  (6  —  a)  region 


(6  —  a) 


P(Xi ) 

[p(xi)dx  =  /V 


—  a) 


dx  =  1 


(2.30) 

(2.31) 


which  is  required  of  p(ay)  (as  a  valid  PDF). 

A  development  of  this  probability  theory  is  contained  in  [5,  section  3.4],  with 
more  detailed  mathematical  proofs  developed  in  [6].  The  definition  concludes  with  a 
measure  of  the  variance  in  I  that  reduces  as  N  increases,  is  detailed  in  Equation  (2.32). 

= *  /  (if _  *) p(x)dx  <2-32) 
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2.5.3  2-Dimensional  Monte  Carlo  Integration.  Extending  the  technique  into 
two  dimensions  is  a  trivial  extension  of  Equation  (2.29), 


I 


f(x,y)dxdy 


(2.33) 


f{xuVi) 
N  P(xi’Vi) 


(2.34) 


2.5.4  Monte  Carlo  Integration  Example.  As  an  example,  let  f(x)  =  5x4  and 
p(x)  =  1  over  the  domain  [0, 1],  which  implies  the  following  integral  and  its  estimate 
from  Equation  (2.29): 


I 

I 


N 


N 


5  a;1 


2=1 


1 


(2.35) 

(2.36) 


The  variance  from  Equation  (2.32)  becomes 

a2  =  Jil{5xl~1)2dx  =  ^-  <2'37) 

Figure  2.16  plots  the  integral  estimate  of  Equation  (2.35)  using  Equation  (2.36)  and 
shows  the  accuracy  of  the  estimate  increasing  with  the  number  of  samples,  N ,  whilst 
remaining  inside  the  bounds  of  the  variance  estimate  calculated  with  Equation  (2.37). 

The  usefulness  of  Monte  Carlo  integration  will  become  evident  in  Section  3.3. 


2.6  Spherical  and  Hemispherical  Geometry 

Figure  2.2  defined  the  geometry  of  the  BRDF  parameters.  Several  concepts 
and  data  display  methods  will  be  used  in  this  thesis,  and  this  section  discusses  the 
utility  of  this  method.  The  BRDF,  after  ignoring  the  wavelength  domain,  is  a  four 
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Figure  2.16:  Monte-Carlo  integration  example.  The  variance 
of  the  solution  decreases  as  the  number  of  samples  increases, 
staying  within  the  variance  bounds. 

dimensional  function.  It  contains  two  domains  (incident  and  reflected  hemispheres) 
of  two  dimensions  each.  Each  hemisphere  has  two  dimensions,  and  is  represented  by 
the  spherical  geometry  variables  <fi  £  [0  :  27t]  and  6  £  [0  :  7r/2]  which  map  out  the 
hemisphere. 

A  solid  angle,  12,  is  defined  by  the  area  subtended  on  the  hemisphere  by  the 
angular  region  of  interest,  and  has  units  of  steradian  (sr).  Assuming  the  hemisphere’s 
radius  is  one,  the  surface  area  of  a  hemisphere  is  then  27rr2  =  27T,  implying  there 
are  47t  steradians  in  a  whole  sphere.  The  area  creating  the  solid  angle  can  be  of  any 
shape,  only  the  magnitude  of  the  area  is  important. 

2.6.1  Hemispherical  Integration.  In  most  cases,  there  is  a  requirement  to  in¬ 
tegrate  the  incident  irradiance,  at  a  point  on  the  model’s  surface,  over  the  hemisphere 
to  calculate  the  reflected  radiance  in  the  direction  of  the  observer.  When  numerically 
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calculating  a  hemispherical  integral,  both  domains  are  treated  as  if  they  are  in  carte¬ 
sian  form.  Thus,  the  differential  solid  angle,  dco,  uses  the  Jacobian  transform  from 
spherical  to  cartesian  coordinates.  When  integrating  the  BRDF  over  the  hemisphere, 
the  differential  solid  angle  (centered  on  the  direction  0), 

dujQ  =  sin  (9)d6d<p  (2.38) 


transforms  the  spherical  integral  into  a  form,  which  the  computer  can  calculate, 


/n  27T 

f(Q)dujQ  =  I 

Q  '  ° 


rn/2 


f(0,  </>)  sin {6)d6d(j). 


(2.39) 


2.7  Chapter  Summary 

The  background  theory  and  concepts  discussed  above  will  be  expanded  and 
applied  to  the  IR  simulation  application  in  Chapter  III.  Chapter  III  also  includes  the 
methods  and  algorithms  used  to  compute  the  IR  signature  simulation. 
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III.  Development  of  the  Infrared  Signature  Simulation 


This  chapter  describes  the  method,  procedures  and  ideas  that  will  be  simulated 
in  Chapter  IV.  The  background  theory  presented  in  Chapter  II  is  expanded  in 
detail  and  is  focussed  towards  the  IR  signature  modelling  application. 

3.1  Infrared  Signature  Components 

An  IR  signature  is  the  signal  that  is  emitted  or  reflected  from  a  target  (aircraft) 
across  the  whole  IR  spectrum.  It  contains  thermal  self  emission  from  the  target  and 
reflections  from  the  background.  The  background  consists  of  earth  and  sky  shine, 
solar  reflections  and  atmospheric  path  radiance;  all  of  these  are  reflected  in  some  part 
from  the  target  towards  the  detector.  The  target’s  composite  signal  is  then  attenuated 
by  the  atmosphere  due  to  scattering  and  absorption. 

Additionally,  the  background  itself  provides  a  direct  component.  The  back¬ 
ground  that  is  in  the  field  of  view  (FOV)  of  the  detector  contributes  through  the  path 
radiance  to  infinity  and  the  path  radiance  between  the  detector  and  target. 

For  single-pixel  detectors,  the  whole  FOV  is  integrated  into  a  single  output,  and 
for  the  target  to  be  detected,  it  must  have  a  greater  magnitude  than  the  background 
(or  possibly  less  in  a  negative-contrast  system).  However,  for  imaging  detectors, 
there  is  spatial  information  which  can  be  used  to  detect  targets.  Adjacent  pixels  are 
compared  to  determine  contrast,  finding  the  edges  of  the  target.  If  the  resolution 
is  sufficient,  then  algorithms  can  classify  the  type  of  target.  However,  this  process 
requires  a  larger  signal-to-noise  ratio  (SNR)  than  pure  detection. 

3.2  Aircraft  Model 

The  aircraft  or  target  being  modelled  must  be  represented  quantitatively  in 
some  manner  to  allow  algorithms  and  mathematical  expressions  to  be  applied  and 
calculated,  respectively.  LCAIR  uses  a  facetised  approach,  where  the  aircraft  surfaces 
are  broken  (facetised  or  tesselated)  into  small  flat  polygons  of  at  least  three  coplanar 
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vertices.  The  number  of  facets  is  arbitrary,  however,  the  more  facets,  the  better 
the  approximation  of  the  real  aircraft.  This  facet-based  approach  is  more  commonly 
known  as  a  wireframe  model. 

3. 2. 1  Wireframe  Model.  Each  component  of  the  real  aircraft  or  target  must 
be  facetised.  The  resulting  model  should  accurately  represent  the  geometry  for  the 
simulation  result  to  be  representative  of  the  real  world  phenomenon.  To  achieve  this, 
each  facet  must  be  continuous  to  the  next  so  there  are  no  gaps  where  sample  rays 
could  miss  an  intersection  between  two  facets.  The  wireframe  facets  are  defined  by 
three  or  four  vertices  (for  a  triangle  or  a  quadrilateral,  respectively)  in  a  single  2D 
(coplanar)  plane.  The  units  are  not  important,  however,  to  aid  in  the  calculation 
of  radiance  later  (with  units  of  watts  per  square  meter  per  steradian ),  meters  are 
chosen  as  the  base  unit.  In  addition  to  the  vertices,  each  facet  must  have  several 
other  parameters  defined.  The  wireframe  file  contains  all  the  information  describing 
the  aircraft  (target)  that  is  required  by  the  propagation  algorithm. 

3. 2. 1.1  Normal.  The  facet  normal  defines  the  direction  that  the  facet 
emits  and  reflects  radiation.  A  condition  of  this  is  that  all  facets  are  only  one-sided. 
The  model  must  be  continuous  and  enclose  the  volume  of  the  aircraft  completely, 
and  all  aircraft  parts  (especially  thin  parts)  must  be  facetised  into  a  volume  with  the 
normals  facing  out. 

3. 2. 1.2  Center.  The  center  of  the  facet  is  calculated  by  the  arithmetic 
mean  of  the  vertices  in  all  three  dimensions.  The  center  is  used  to  calculate  the 
spherical  angles  (</>  and  0)  to  and  from  the  other  facets.  The  irradiance  and  exitance 
distribution  is  assumed  constant  across  an  individual  facet,  with  the  center  being  used 
for  the  geometric  point  of  calculation. 

3.2. 1.3  Reference  Vector.  The  azimuthal  domain  is  expressed  from 
—  180°  to  +180°;  the  reference  vector  divides  the  positive  and  negative  angles.  The 
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normal 


vector 


positive  vector 


Figure  3.1:  Sketch  of  the  definition  of  the  positive  and  neg¬ 

ative  azimuth  domain  values.  The  positive  half  always  has  the 
positive  vector,  as  this  defines  a  right  handed  coordinate  system. 


reference  vector  is  calculated  from  the  facet  center  to  the  second  vertex,  and  is  nor¬ 
malized  to  give  a  unit  vector. 


3.2. 1-4  Positive  Vector.  The  positive  vector  defines  the  positive  180° 
of  azimuth  and  lies  along  the  [(f)  =  90°,  6  =  90°]  vector.  It  is  calculated  by  a  cross 
product  to  create  a  right  handed  coordinate  system  with  the  reference  vector  and  facet 
normal.  Figure  3.1  defines  the  positive  and  negative  azimuthal  values  in  relation  to 
the  reference  and  positive  vectors. 


3.2. 1.5  Facet  Type.  Each  facet  is  coded  to  indicate  what  it  represents, 
whether  a  model  (f),  sky  (s)  or  terrain  (t)  facet. 

3. 2. 1.6  Material  type.  Surface  treatment  of  the  facet  is  very  important 
as  it  affects  the  reflectance  distribution  of  the  facet.  The  material  type  is  coded  as 
the  propagation  algorithm  must  know  what  the  surface  treatment  is.  Material  types 
available  for  simulation  are  gloss  black  painted  aluminum,  bare  aluminum,  and  a  per- 
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fectly  diffuse  surface.  Other  surfaces  which  are  available,  although  not  simulated  are 
FS3411  Flat  Blue  Camouflage  and  FS3464  Gloss  White  as  the  Sandford- Robertson 
model  parameterizations  are  listed  in  [23].  The  Federal  Standard  Paint  scheme  (ab¬ 
breviated  as  FSxxxx)  is  an  attempt  to  standardize  colors  in  the  United  States  Federal 
Government  to  provide  a  way  for  engineers  and  manufacturers  to  reference,  design 
and  manufacture  to  one  common  paint  scheme.  FED-STD-595,  [1],  is  the  controlling 
document  for  the  standard. 

3.2. 1.7  Temperature.  As  discussed  in  Chapter  II,  the  temperature 
of  a  surface  affects  the  blackbody  radiation  of  a  surface.  Each  facet  must  have  a 
temperature  assigned  in  Kelvin. 

3.2.2  Blender  3D  Program.  The  user  needs  a  tool  to  create,  edit  and 

manipulate  the  wireframe  model  with  ease.  The  software  program  Blender3D  was 
used  to  create  the  wireframe  and  assign  the  facet  type,  material  type  and  temperature 
to  each  facet. 

An  off-the-shelf  product  was  primarily  chosen  for  ease  of  use,  but  also  because 
it  comes  with  an  import / export  utility  that  can  import  many  types  of  3D  wireframes 
that  are  available  for  aircraft.  Autodesk  3D  Studio  Max@  hies  are  easily  imported. 
Additionally,  once  the  wireframe  model  is  loaded  inside  Blender3D,  the  export  utility 
is  used  to  produce  a  data  hie  for  the  Matlab®algorithm  to  read,  containing  the  vertex 
and  label  format  information  for  each  facet.  The  data  hie  in  Blender3D  is  exported 
as  a  Wavefront  (.obj)  hie.  The  settings  required  for  this  are  shown  in  Figure  3.2. 

Each  facet  is  coded  with  the  parameters  that  were  dehned  in  Section  3.2.1. 
The  held  options  are  listed  in  Table  3.1,  and  the  label  must  conform  to  the  format 
in  Figure  3.3.  An  example  is  included  for  a  bare  aluminum  model  facet  with  a 
temperature  of  300K.  Figure  3.4  contains  a  screenshot  of  Blender3D,  which  shows 
the  facet  label  format  entry  position  in  the  top  left  of  the  textbox. 
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Mesh  Options... 

Apply  Modifiers 

Triangulate 
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High  Quality  Normals 

UVs 

Materials 
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Context... 

Selection  Only 

All  Scenes 

Animation 

Copy  Images 

Grouping... 

Objects 

Object  Groups 

Material  Groups 

Figure  3.2: 


Screenshot  of  the  Wavefront  (.obj)  export  settings  in  Blender3D. 


■{facet  type}* {material  type}* {temperature} *{Blender3D  misc} 
f *bareAl*300* 

Figure  3.3:  The  facet  naming  format  must  be  typed  exactly  so  the  Matlab®  import 
function  can  read  the  applicable  fields.  This  example  is  for  a  bare  aluminum  model 
facet  at  a  temperature  of  300K. 


Figure  3.4:  Screenshot  of  the  cylinder  model  in  Blender3D.  The  program  makes 

editing  and  manipulating  the  wireframe  model  very  easy. 
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Table  3.1:  Facet  label  attributes  for  facet  type, 

material  type  and  temperature.  These  must  be  en¬ 
tered  exactly  into  the  textbox  in  Blender3D,  or  the 
Mat  lab®  import  function  will  crash. 


Label  Field 

Format  Options 

Detail 

Facet  Type 

f 

model  facet 

s 

sky  facet 

t 

terrain  facet 

Material  Type 

glossBk 

Gloss  black  painted  aluminum 

bareAl 

Bare  aluminum 

diffuse  1.0 

Perfectly  diffuse 

(Lambertian)  blackbody,  e  =  1.0 

diffuse0.5 

Perfectly  diffuse  greybody,  e  =  0.5 

diffuseO.O 

Perfectly  diffuse  reflector,  e  =  0.0 

modt^ 

MODTRAN  generated  path  radiance 

Temperature 

t  [k; 

Temperature  of  facet  in  Kelvin 

Blender3D  misc 

{blank} 

Blender3D  adds  information  to  the  label. 
All  information  to  the  right  of  the 
third  asterisk  (*)  is  ignored 
by  the  Matlab®import  function. 

3.2.3  Simulation  Models.  Three  models  are  simulated  in  Chapter  IV;  the 
first  is  a  cavity  blackbody  source.  The  second  is  a  cylinder  with  wings,  representing 
a  notional  aircraft  shape.  It  is  simulated  in  an  IR  range  environment  to  demonstrate 
the  applicability  and  accuracy  of  the  various  algorithms.  The  final  model  is  based  on 
an  Airbus  A340-300,  and  is  simulated  at  a  generic  airport  to  calculate  the  apparent 
irradiance  at  a  detector  located  outside  the  airport  boundary.  Other  models  are 
included  in  Appendix  A,  where  analytical  test  cases  are  performed  to  quantitatively 
validate  the  algorithms. 

3.2.3. 1  Artificial  Cavity  Blackbody  Source.  Figure  3.5  shows  a  cavity 
model  displayed  in  Matlab®.  The  cavity  model  is  facetised  into  400  facets  and  has 
a  radius  of  0.5m.  During  the  simulation,  seven  other  cavities  were  simulated  with 
varying  aperture  sizes.  The  eight  different  aperture  sizes  are  created  by  removing 
more  and  more  facets  from  the  wireframe  model.  The  number  of  facets  removed  are; 


Figure  3.5:  The  artificial  cavity  blackbody  source  is  modelled 
as  a  facetised  sphere  with  400  facets.  In  this  case,  two  facets  are 
removed  to  create  the  aperture. 

2,  4,  8,  16,  24,  36,  60,  and  80.  Figure  B.l  shows  the  different  cavity  aperture  sizes. 
The  cavity  is  set  to  800  Kelvin  and  has  an  internal  emissivity  of  0.5. 

3. 2. 3. 2  Cylinder  Model.  The  cylinder  model  was  built  in  2005  by 
Capt  Bortle  as  a  physical  model  to  validate  the  original  LCAIR  code  written  by 
Capt  Martinez  in  2005.  The  model  is  constructed  of  aluminium  and  has  dimensions 
as  detailed  in  Figure  3.6.  Detailed  construction  drawings  can  be  found  in  [3]  and 
photographs  of  the  model  are  included  in  Figures  3.7  and  3.8.  Figure  3.9  shows  the 
cylinder  model  in  the  Blender3D  program.  The  cylinder  model  is  simulated  in  two 
configurations,  bare  aluminum  and  painted  with  gloss  black  paint.  Table  3.2  lists  the 
temperatures  that  were  measured  in  the  IR  range  simulations  by  Bortle. 

3. 2. 3. 3  Airbus  A340-300.  The  third  model  to  be  simulated  is  the 
Airbus  A340-300.  This  aircraft  was  chosen  as  a  typical  large  commercial  airliner. 
A  sample  wireframe  model  was  imported  into  Blender3D.  From  there,  the  A340-300 
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Table  3.2:  Cylinder  Model  Surface  Temperatures. 

Reproduced  from  [3]. 


Position 

Gloss  Black 

Temperature  [K] 

Bare  Aluminum 

Temperature  [K] 

Nose 

292.1 

301.7 

Front  Fuselage 

293.1 

302.1 

Rear  Fuselage 

294.2 

303.2 

Left  Wing  Front 

301.6 

303.6 

Left  Wing  Rear 

298.2 

303.6 

Right  Wing  Front 

304.2 

303.7 

Right  Wing  Rear 

299.1 

303.7 

$3” 


-»< - ><- 


6”  3”  6” 

Front  View 


✓ 

3” 

1 

|o  3” 

V 

✓ 

3”  < - 

— > 

< — 

— > 

iN1 


16” 
Side  View 


3” 


16” 


Figure  3.6:  3-View  drawing  of  the  cylinder  model,  showing  the  dimensions  from 

the  front,  side  and  top  views.  Adapted  from  [3,  Figure  3.1]. 
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(a)  (b) 


Figure  3.7:  A  photo  of  the  cylinder  model  showing  (a)  the  internal  construction  of 
the  cooling  coils  and  the  thermal  electric  heater  input  wires,  and  (b)  another  photo 
showing  the  cylinder  model  with  the  cover  attached.  In  this  figure  the  model  is  not 
painted,  which  is  the  surface  treatment  for  the  bare  aluminum  simulations  in  the  IR 
Range.  Reproduced  from  [3,  Figure  3.2]. 


Figure  3.8:  Photo  of  the  cylinder  model  with  the  gloss  black  paint  surface  treatment 
applied.  Note  the  underside  of  one  wing  is  left  bare  to  provide  some  contrast  in  the 
IR  range  simulations.  Reproduced  from  [3,  Figure  3.12], 
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Figure  3.9:  A  screenshot  of  the  cylinder  model  in  the 

Blender3D  program,  viewed  from  different  orientations. 


model  is  treated  exactly  the  same  as  the  cylinder  model  with  respect  to  the  algorithms 
and  simulation  procedures.  Figure  3.10  shows  the  aircraft  model  in  Matlab®. 

During  the  simulations,  the  aircraft’s  surface  treatment  will  be  bare  aluminum, 
to  approximate  the  specularity  of  commercial  aircraft  surface  finishes.  The  aircraft 
fuselage,  wings,  stabilizers,  engine  pylons  and  cowls  are  set  to  300K,  the  ambient 
temperature.  The  majority  of  the  thermal  signature  is  simulated  to  be  emitted  from 
the  engine  exhaust  plume  and  the  visible  hot  parts  in  the  turbine  cavity.  The  exhaust 
plume  is  modelled  as  a  simple  cone  at  a  temperature  of  500K.  Figure  3.11  shows  the 
exhaust  temperature  profile  for  a  turbofan  engine,  and  is  assumed  to  be  similar  to 
the  A340-300  engine.  From  the  figure,  the  500K  profile  would  extend  to  about  20  feet 
from  the  tailpipe  exit  plane;  hence,  the  exhaust  cone  in  the  model  is  set  to  20  feet  in 
length. 

Additionally,  the  engine  hot  parts  are  modelled  by  an  annulus  that  connects  the 
exhaust  cone  to  the  engine  cowl.  This  simulates  the  view  of  the  engine  from  behind 
where  the  observer  can  see  into  the  engine  cavity;  the  temperature  is  set  to  800K. 
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Figure  3.10:  The  Airbus  A340-300  model  is  imported  into  Matlab®and  displayed 
from  three  angles.  Note  the  varying  magnitude  of  the  visual  cross-section  from  the 
three  different  angles. 
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Figure  3.11:  Engine  exhaust  temperature  contours  for  the 

Boeing  707  turbofan  and  turbojet  engines.  Reproduced 
from  [26,  Figure  2-72], 

This  annulus  is  not  visible  forward  of  the  tailpipe  exit  plane;  which  again,  simulates 
the  view  into  the  engine  cavity. 

3.2.4  Hemispherical  View.  The  hemispherical  view  of  any  point  on  the 
model’s  surface  refers  to  the  geometry  in  the  hemisphere  that  is  visible  from  that 
point.  Viewable  objects  include  other  model  surfaces  or  background,  whether  terrain 
or  sky.  The  hemispherical  view  has  information  of  where  (as  a  function  of  </>  and 
9)  and  how  large  (size  of  solid  angle)  each  object  is.  The  view  is  displayed  in  a 
2D  figure  with  the  colormap  referenced  to  the  3D  model  facet  index.  Figure  3.12 
shows  the  hemispherical  view  of  facet  A  from  the  example  geometry  of  Figure  2.2, 
and  Figure  3.13  shows  the  view  of  the  more  complex  geometrical  setup  of  the  cylinder 
model  from  an  arbitrary  facet  on  the  wing  leading  edge. 

Displaying  BRDF  values  over  the  hemisphere  requires  a  similar  2D  figure.  In 
this  case,  one  of  the  two  domains  must  be  fixed,  to  allow  the  function  to  be  evaluated 
across  the  other.  The  most  common  display  is  when  the  reflected  domain  is  fixed 
at  some  observation  position  and  the  BRDF  is  calculated  across  the  whole  incident 


44 


domain;  an  example  with  9i  =  30°,  d%  =  100°  is  shown  in  Figure  3.14.  Recall  that 
under  the  theory  of  reciprocity  from  Section  2. 3. 1.1,  it  does  not  matter  which  domain 
is  fixed;  the  result  will  be  the  same  as  described  in  Equation  (2.6). 

3.3  Integrating  the  BRDF 

Section  2.3.2  discussed  the  different  integral  forms  of  the  BRDF;  the  most  useful 
one  is  the  DHR,  where  the  BRDF  is  integrated  over  all  exitant  angles.  This  hemi¬ 
spherical  integration  will  be  required  by  the  propagation  algorithms  of  Section  3.5. 
Hemispherical  integration  is  required  when  integrating  the  incident  irradiance  of  any 
facet.  The  BRDF  must  also  be  hemispherically  integrated  to  scale  that  incident  ir¬ 
radiance.  Computer  simulations  require  all  calculations  (including  integrals)  to  be 
computed  numerically  (discretely).  The  caveat  to  this  approach  is  that  the  function 
must  be  suitably  sampled  during  the  integral  to  avoid  aliasing  problems  that  give  an 
incorrect  answer.  This  sampling  issue  is  the  same  as  discussed  in  Section  2.5,  where 
Figure  2.16  shows  the  aliasing  error  decreasing  as  the  number  of  samples  increases. 

The  BRDF  function  must  be  suitably  sampled  during  integration  to  ensure  that 
the  function’s  characteristics  are  observed  throughout  the  domains.  The  Sandford- 
Robertson  BRDF  model,  being  a  fractional  specularity  model,  assumes  that  the  sur¬ 
face  exhibits  some  specular  behavior.  It  is  this  specular  nature  that  requires  careful 
sampling  when  integrating  the  function.  Figure  2.9  demonstrated  the  specular  be¬ 
havior  of  the  surface  where  the  diffuse  floor’s  magnitude  is  10-2  and  the  specular 
peak  is  10+4.  This  is  a  very  large  range  across  the  function’s  domain.  The  region 
where  the  large  values  (10+1  and  above)  exist  is  also  over  a  small  angular  region  of 
approximately  5°  in  9r,  which  demonstrates  the  specularity  of  the  aluminum  surface. 
Note  the  gloss  black  surface,  although  with  a  considerably  lower  total  reflectance  than 
aluminum,  is  also  highly  specular;  the  angular  region  of  large  values  (10°  and  above) 
is  also  5°  in  9r. 

Having  to  numerically  integrate  this  function  creates  a  computational  problem. 
A  design  goal  of  any  simulation  is  low  computation  time,  hence,  the  number  of  samples 
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(a) 


Azimuth  (degrees) 

(b) 


Figure  3.12:  The  geometry  of  a  simple  setup  is  shown  in  (a).  The  colormap  refer¬ 
ences  to  the  facet  index  number.  Note  the  region  in  <fr  and  6  that  each  of  the  other 
facets  occupy  in  (b).  Also,  observe  the  spherical-to-cartesian  conversion  effect,  which 
makes  square  facets  have  curved  edges  in  the  hemispherical  view  representation. 
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(a) 


Azimuth  <|)  (degrees) 

(b) 

Figure  3.13:  (a)  Geometry  of  the  cylinder  model.  The  colormap  is  referenced  to  the 

facet  index  number,  (b)  The  hemispherical  view  of  a  wing  leading  edge  facet  (#1154) 
from  the  cylinder  model.  Note  the  region  where  the  model  facets  exist,  everywhere 
else  is  background. 
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Figure  3.14:  (a)  The  BRDF  is  represented  in  2D  format.  This  data  structure  is  a 

2D  matrix  that  contains  values  as  a  combination  of  azimuth  and  elevation  dimensions, 
(b)  The  BRDF  is  shown  in  3D  here  for  clarity.  Note  the  log  scale  in  the  y-axis  and 
colormap. 
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Figure  3.15:  Integrated  BRDF  of  bare  aluminum 

(6.j  =  30°,  fa  =  100°)  in  the  long  wave  IR  (8  —  12/im)  with 
varying  sample  resolution  in  the  azimuth  and  elevation  do¬ 
mains.  The  DHR  (reflectance)  for  LWIR  bare  aluminum  with 
6i  =  30°  is  0.93. 

per  integration  domain  must  be  kept  as  low  as  possible  to  reduce  computation  time. 
However,  this  need  is  balanced  against  the  requirement  to  adequately  sample  the 
function.  Recall  from  Equation  (2.8),  the  integral  of  the  BRDF  over  all  reflected 
angles  must  equal  the  reflectance  at  that  reflected  elevation  angle.  Simulations  were 
performed  by  calculating  the  DHR  from  the  BRDF  with  varying  resolution  across 
the  azimuth  and  elevation  domains  to  determine  the  sampling  resolution  required  to 
achieve  an  accurate  result,  whilst  keeping  the  computation  time  minimized. 

Figure  3.15  shows  the  DHR  result  for  sampling  resolutions  of  100  to  700  samples 
per  dimension.  Note  that  the  number  of  samples  is  per  dimension,  thus,  to  integrate 
the  BRDF  across  or  and  6r  (two  dimensions)  requires  N2  samples.  From  Figure  3.15, 
500  samples  per  dimension  are  required  for  an  accurate  result,  thus,  5002  =  250,  000 
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samples  are  actually  computed.  This  number  of  samples  is  prohibitive  for  a  timely 
simulation,  even  with  a  single  bounce  simulation,  as  each  facet’s  observed  radiance 
takes  minutes  to  calculate.  The  relatively  simple  cylinder  model  has  around  2,000 
facets;  with  several  minutes  per  facet,  a  single  bounce  simulation  will  take  days! 
Thus,  an  optimized  result  must  be  found. 

As  an  aside,  notice  in  the  result  of  the  BRDF  integrated  DHR  measurement 
in  Figure  3.15;  the  solution  converges  on  0.917,  which  is  lower  than  the  empirically 
measured  DHR  of  0.93  from  Figure  2.7.  This  difference  results  from  the  fitment  of 
the  BRDF  model  parameters  to  the  empirical  BRDF  curves.  It  can  be  seen  that 
the  parameterized  BRDF  curves  in  Figure  2.9  fall  under  the  empirical  measurement 
everywhere  except  the  specular  peak.  When  the  parameterized  curve  is  integrated, 
the  result  is  less  (owing  to  the  lesser  area  under  the  curve).  The  opposite  can  be 
demonstrated  for  the  gloss  black  case  (Figure  2.12),  where  a  significant  proportion 
of  the  specular  peak  remains  above  the  empirical  curve,  whilst  the  diffuse  floors  are 
relatively  equal. 

Importance-sampled  Monte  Carlo  integration  and  numerical  quadrature  will 
now  be  investigated  to  search  for  the  optimum  numerical  integration  technique  for 
this  IR  signature  simulation  application. 

3.3.1  Importance  Sampled  Integration.  The  blind  numerical  integration 
technique  simulated  above  has  no  a  priori  information  about  the  BRDF  distribution. 
In  this  case,  the  majority  of  samples  are  wasted,  especially  when  the  BRDF  resembles 
a  delta  function  where  the  values  of  any  appreciable  magnitude  are  so  closely  located. 

A  more  efficient  calculation  method  is  to  sample  the  function  with  a  density 
that  is  proportional  to  the  value  of  the  function,  where  the  peaks  are  sampled  many 
more  times  than  the  lower  values.  This  method  places  importance  on  the  larger 
values  of  the  function  as  they  have  a  greater  bearing  on  the  result.  Recall  that  when 
numerically  integrating,  a  weighting  factor  l/p(aq)  must  be  applied  to  the  function’s 
sampled  value.  The  value  of  p(xi )  is  the  probability  that  the  function  will  be  sampled 
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Figure  3.16:  Integration  with  importance  sampled  domain. 

The  p(xi)  function  is  representative  only,  and  it  is  not  meant 
to  specifically  depict  a  piecewise  function.  The  choice  of  the 
distribution  for  p(xi)  is  an  important  component  of  achieving 
the  fastest  convergence. 

at  that  point,  and  the  higher- probability  sample  values  must  be  weighted  by  a  larger 
factor  to  account  for  them  being  closely  spaced.  The  region  of  the  function  with  a 
higher  probability  of  being  sampled  will  then  have  more  samples  compared  to  the 
other  regions,  as  shown  in  Figure  3.16,  where  the  error  introduced  by  the  rectangular 
approximation  is  much  less. 

The  choice  of  distribution  for  p(xi)  should,  optimally,  be  related  to  the  function 
itself.  The  perfect  p(xi )  is  then, 


p(Xi) 


1/0*01 

f  f(x)dx 

D 


(3.1) 


where  the  function  itself  (when  converted  to  a  probability  distribution  function)  in¬ 
tegrates  to  one.  Whilst  this  is  the  perfect  p(xi),  which  will  result  in  the  variance  of 
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Equation  (2.32)  equaling  zero,  it  also  requires  knowledge  of  the  integral  of  the  BRDF, 


J  f(x)dx  =  /  (3.2) 

D 

which  is  usually  the  object  of  the  numerical  integration  estimate,  I. 

3.3.2  Monte  Carlo  Integration.  Monte  Carlo  integration,  introduced  in  Sec¬ 
tion  2.5,  is  the  term  for  the  technique  where  the  samples  along  one  or  more  dimensions 
are  generated  with  a  certain  distribution,  p(xi). 

From  the  previous  section,  the  exact  integral,  I,  is  achieved  with  the  perfect 
distribution  from  Equation  (3.1).  However,  the  key  is  being  able  to  generate  samples 
from  the  p(aq)  distribution.  The  most  efficient  way  is  to  use  the  inverse  Cumulative 
Distribution  Function  (CDF)  method,  where  a  uniform  random  variable  is  trans¬ 
formed  into  the  distribution  of  choice.  This  method,  however,  requires  an  expression 
for  the  CDF  of  the  random  variable,  which  is  based  on  the  BRDF  of  the  specific 
surface  treatment.  Section  2.3.4  detailed  the  method  of  calculating  the  BRDF  value 
using  the  four  parameters  and  the  equations  defined  in  [23],  however,  the  complexity 
of  the  BRDF  model  prevented  an  analytical  expression  of  the  CDF  (matching  the 
reflectance  distribution  of  the  BRDF)  from  being  developed.  In  theory,  the  analytical 
CDF  is 


CDF(9i,<j)i,Qr,(f>r ) 


|  fr  (ftii  0i>  9r,  (f>r)  | 

§  fr  (C; ,  (f>i ,  0T ,  (pr  )  (10f  d(f)r 


dOrd(f)r. 


(3.3) 


\D  / 

Other  simpler  BRDF  models  are  able  to  be  integrated  analytically;  however,  the 
Sandford-Robertson  model  is  being  used  in  this  case. 


An  analytical  expression  for  the  CDF  is  required  to  generate  samples  quickly  and 
efficiently.  A  method  of  developing  samples  to  a  distribution  without  an  analytical 
CDF  is  called  rejection  sampling.  Rejection  sampling  is  where  uniform  samples  are 
generated  in  N+l  domains,  where  N  is  the  number  of  domains  in  p(x).  The  samples 
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are  then  tested  against  the  required  distribution,  where  if  they  fall  under  the  desired 
p(x)  curve,  they  are  accepted.  A  more  detailed  description  of  this  technique  and  other 
generation  techniques  can  be  found  in  [6].  Despite  this  technique  being  available  to 
generate  samples,  the  delta  function-like  distribution  of  values  in  the  BRDF  makes 
this  method  extremely  inefficient,  and  it  is  not  chosen.  The  location  of  the  BRDF 
peak  is  different  for  each  evaluation  of  the  reflected  radiance  integral,  which  would 
make  even  a  precalculated  look-up  table  method  inaccurate. 

3.3.2. 1  Importance  Sampled  Normal  Distribution.  Although  a  BRDF- 
based  CDF  was  not  found,  an  attempt  to  integrate  using  a  normally  distributed 
approximation  was  performed.  The  normal  distribution  was  chosen  as  its  CDF  is  well 
known;  thus,  samples  can  be  calculated  easily.  Figure  3.17  shows  the  convergence  of 
/  using  the  normal  distribution  approximation.  The  integral  approximation  begins  to 
converge;  however,  at  250,000  samples,  the  technique  has  still  not  achieved  a  suitable 
accuracy.  This  number  of  samples  is  clearly  time  prohibitive,  and  this  technique  is 
rejected. 


3. 3. 2. 2  Stratified  Uniform  Distribution.  Another  adaptation  to  the 
Monte  Carlo  case  is  stratified  sampling,  where  the  function’s  domain  is  divided  into 
several  regions  where  different  sampling  resolutions  can  be  used  in  each.  The  regions 
can  be  sampled  with  uniform  distributions  as  depicted  in  Figure  3.16. 

The  first  region  is  defined  as  the  function’s  peak  centered  on  the  specularly 
reflected  point  with  a  width  of  20°  in  elevation  and  40°  in  azimuth.  However,  the 
statistical  nature  of  a  uniform  random  variable  causes  the  technique  to  break  down 
as  the  uniformly  generated  variable  is  not  purely  uniform.  The  effect  of  this  ap¬ 
proximation  and  lack  of  uniformity  is  demonstrated  in  Figure  3.18,  where  a  suitable 
accuracy  is  not  reached  after  500,000  samples.  Special  algorithms  exist  for  generat¬ 
ing  better  distributed  uniform  random  variables,  however,  they  were  not  investigated 
as  the  numerical  quadrature  section  below  proves  to  be  far  superior.  Examples  of 
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Figure  3.17:  BRDF  integration  with  a  normal  distribution 

approximation.  Note  how  after  250,000  samples,  the  approxi¬ 
mation  has  not  suitably  converged  to  the  correct  result. 

these  optimized  distributions  can  be  found  in  [6] ,  where  one  technique  is  based  on  the 
Halton  sequence  of  pseudo  random  numbers. 

3.3.3  Stratified  Numerical  Quadrature.  The  advantage  of  statistical  tech¬ 
niques  comes  from  the  law  of  large  numbers,  and  is  only  realized  when  the  number  of 
samples  becomes  large.  It  has  been  demonstrated  that  even  500,000  samples  is  not 
suitably  large  enough  when  using  this  technique  for  this  application. 

A  variation  of  the  stratified  uniform  distribution  is  simply,  stratified  numerical 
quadrature.  In  this  case,  the  same  regions  as  Section  3. 3. 2. 2  were  chosen  and  inte¬ 
grated  with  numerical  quadrature.  Figure  3.19  shows  the  convergence  of  the  peak 
region  integral  approximation,  whilst  Figure  3.20  shows  the  convergence  of  the  outer 
region  integral. 
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Figure  3.18:  BRDF  integration  with  a  stratified  uniformly 

distributed  approximation.  Note  how  after  500,000  samples,  the 
approximation  has  not  suitably  converged  to  the  correct  result. 
A  single  dot  represents  one  simulation  at  that  number  of  sam¬ 
ples.  If  the  dots  are  widely  spaced,  then  subsequent  simulations 
provide  different  results,  and  hence,  a  large  variance  in  the  so¬ 
lution. 


This  method  is  far  superior  to  those  previously  simulated.  The  advantage  of 
stratified  numerical  quadrature  is  that  the  peak  is  sampled  at  a  higher  resolution  than 
the  rest  of  the  domain.  To  sample  at  this  resolution  across  the  whole  domain  would 
require  many  samples,  which  was  shown  in  Section  3.3. 


3.3.4  Geometry  Sampling.  In  addition  to  integrating  the  BRDF  correctly, 
the  hemisphere  must  also  be  sampled  at  a  high  enough  rate  to  locate  the  edges  of  each 
facet.  Figure  3.21  shows  the  number  of  samples  required  to  sample  the  facet’s  radiance 
when  integrating  across  the  BRDF.  The  figure  shows  that  around  200  samples  per 
dimension  are  required  for  a  correct  result.  This  sampling  rate  is  required  as  each 
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Figure  3.19:  Convergence  of  the  peak  region  (40°  x  20°)  is  achieved  much  faster 

with  numerical  quadrature.  After  only  approximately  80  samples  per  dimension,  the 
correct  result  is  achieved. 


Figure  3.20:  Convergence  of  the  outer  region  is  achieved  much  faster  with  numerical 
quadrature.  After  approximately  50  samples  per  dimension,  the  correct  result  is 
achieved;  note  the  resolution  in  the  y-axis  scale. 
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Figure  3.21:  Additional  to  the  BRDF,  the  hemispherical  ge¬ 
ometry  must  be  sampled  to  accurately  calculate  the  solid  angles 
that  each  facet  occupies. 


sample  approximates  the  radiance  in  that  solid  angle.  The  more  samples,  the  smaller 
the  solid  angles  become,  increasing  the  accuracy  of  the  estimated  radiance. 

3.4  Background  Construction 

In  addition  to  the  target’s  self  emission,  there  is  also  a  significant  contribution 
to  the  signature  from  the  background,  especially  when  a  highly  reflective  surface 
treatment  is  applied.  The  sources  of  this  background  component  include  terrain 
reflection  and  self  emission,  sun  reflection  and  atmospheric  path  radiance.  For  this 
reason,  a  detailed  background  model  is  required  for  the  simulation  to  be  accurate, 
even  for  laboratory  simulations. 

3-4-1  IR  Range.  The  IR  Range  at  the  Sensors  Directorate,  AFRL,  is  a  tem¬ 
perature  and  atmosphere  controlled  enclosure  designed  to  simulate  conditions  up  to 
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Figure  3.22:  AFRL’s  IR  Range  is  developed  in  Blender3D  and 
imported  into  Matlab®.  The  cylinder  model  is  included  in  the 
figure  for  a  size  comparison. 

35,000  feet  in  altitude.  Each  wall  is  cooled  with  liquid  nitrogen  and  the  upper  half  can 
be  controlled  separately  to  the  lower  half;  this  allows  a  sky  and  earth  to  be  simulated. 
Cooling  with  liquid  nitrogen  provides  a  temperature  controlled  environment  and  also 
minimizes  the  thermal  emission  of  the  background. 

To  represent  the  IR  Range  background  in  LCAIR,  a  cube  with  sides  of  length, 
3m,  is  facetised  (tesselated)  into  160  facets,  as  shown  in  Figure  3.22,  where  the  cylinder 
model  is  also  included  for  scale  comparison.  Each  facet’s  material  property  is  set  to 
diffuse  emission  and  emissivity  is  set  to  96%.  From  data  included  in  [3],  the  upper 
walls  of  the  range  are  set  to  263K  and  253K  to  simulate  sky  and  terrain,  respectively. 


3-4-2  World  Sphere.  For  a  simulation  to  be  accurate,  it  must  be  represen¬ 
tative  of  the  physical  environment.  The  IR  Range  background  discussed  above  is  a 
suitable  background  for  laboratory  experiments,  however,  it  is  not  as  suitable  for  a 
real  world  simulation.  The  real  world  contains,  earth,  sky,  sun  and  atmosphere  that 
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Figure  3.23:  The  world  sphere  background  model  is  developed 
in  Blender3D  and  imported  into  Mat  lab®.  The  A340  aircraft  is 
included  in  the  figure  to  demonstrate  the  simulation  setup. 

is  contributing  to  the  background  irradiance  onto  the  model.  The  world  sphere  is  de¬ 
signed  to  simulate  sky  and  terrain  when  the  algorithm  does  not  encounter  an  aircraft 
facet.  Figure  3.23  shows  the  world  sphere  with  A340-300  aircraft  inside. 

For  the  A340-300  simulations,  the  ambient  weather  conditions  are  set  to  be 
standard  overcast  day,  at  the  Dallas  Airport  in  Texas,  with  an  air  temperature  of 
27°C  (300K).  An  inland  airfield  was  chosen  to  remove  the  effects  of  water,  and  provide 
solid  ground  as  the  terrain  model. 

3-4-2. 1  Terrain.  The  terrain  is  modelled  as  a  blackbody.  This  deci¬ 
sion  is  based  upon  spectral  radiance  figures  listed  in  the  The  Infrared  Handbook  [26] . 
The  terrain  composition  around  the  airport  is  assumed  to  be  an  even  mix  of  grass 
and  concrete/asphalt  covered  areas.  This  assumption  is  based  upon  a  heavily  popu¬ 
lated  or  industrialized  area  surrounding  the  airport,  which  is  common  of  many  major 
international  airports. 
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Figure  3.24:  Day  and  night  radiances  of  grass-covered  field, 

reproduced  from  [26,  Figure  3.27].  Note  that  the  measured  ra¬ 
diance  matches  the  radiance  of  a  perfect  blackbody  reasonably 
well,  inferring  a  spectral  emissivity  very  close  to  one. 

Figure  3.24  is  reproduced  from  [26]  and  shows  the  spectral  radiance  of  a  grass 
covered  held.  The  blackbody  curves  for  the  ambient  temperatures  are  provided  and 
match  reasonably  well  to  the  measured  radiance;  hence,  grass  will  be  defined  as  a 
perfect  blackbody  at  the  ambient  temperature. 

Concrete  and  asphalt  can  also  be  approximated  by  a  perfect  blackbody,  how¬ 
ever,  the  apparent  temperature  is  higher  than  ambient.  Figures  3.25  and  3.26  show  the 
spectral  apparent  temperature  with  and  without  cloud  cover  for  concrete  and  asphalt, 
respectively.  Spectrally,  the  emission  is  hat,  which  usually  implies  a  greybody;  how¬ 
ever,  no  emissivity  information  is  given.  The  spectral  peaks  between  3  —  4 fim  on  both 
hgures,  with  no  cloud  cover,  are  due  to  scattered  solar  radiation;  however,  it  is  ignored 
in  the  greybody  assumption,  and  an  averaged  brightness  (apparent)  temperature  is 
calculated.  Note  the  ratio  between  ambient  temperature  and  apparent  temperature 
changes  with  and  without  cloud  cover.  Concrete  is  10  Kelvin  higher  than  ambient 
temperature  with  no  cloud  cover  where  asphalt  is  approximately  20  Kelvin  higher. 
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Figure  3.25:  (a)  Spectral  apparent  temperature  for  concrete  with  no  cloud  cover. 

Ambient  temperature  is  300K.  Reproduced  from  [26,  Figure  3-103(a)].  (b)  Spectral 
apparent  temperature  for  concrete  with  full  cloud  cover.  Ambient  temperature  is 
295K.  Reproduced  from  [26,  Figure  3-103(b)]. 
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Figure  3.26:  (a)  Spectral  apparent  temperature  for  asphalt  with  3/10  cloud  cover. 

Ambient  temperature  is  283K.  Reproduced  from  [26,  Figure  3-105(a)].  (b)  Spectral 
apparent  temperature  for  asphalt  with  full  cloud  cover.  Ambient  temperature  is  285K. 
Reproduced  from  [26,  Figure  3- 105(b)]. 
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Under  full  cloud,  concrete  matches  the  ambient  temperature  and  asphalt  is  still  10 
Kelvin  higher. 

Under  the  assumption  of  an  even  mix  between  grass,  and  concrete  and  asphalt 
combined,  with  full  sun,  the  apparent  terrain  temperature  will  be  7.5  Kelvin  hotter 
than  ambient.  Under  full  cloud,  the  apparent  terrain  temperature  will  be  2.5  Kelvins 
hotter.  The  lower  hemisphere  of  the  world  sphere  represents  terrain  and  all  facets  will 
be  set  to  2.5  Kelvin  above  ambient,  302.5  Kelvin. 

3. 4 .2. 2  Sky.  The  sky  is  very  rarely  modelled  as  a  grey-  or  black- 
body.  Atmospheric  transmission  and  radiance  simulations  are  required  to  produce 
the  spectral  radiance  curves. 

The  upper  hemisphere  of  the  world  sphere  represents  the  sky.  The  hemisphere 
is  divided  in  to  10  horizontal  rings  for  elevation  dependence,  which  are  divided  into 
20  facets  themselves,  to  provide  azimuth  dependence  (see  Figure  3.23). 

The  weather  conditions  that  will  be  used  for  the  simulation  are  a  overcast  day 
with  an  ambient  temperature  of  300K.  The  overcast  day  is  chosen  to  remove  the  effect 
of  the  sun.  Including  the  sun  in  the  background  will  influence  the  intensity  of  the 
simulation  by  making  the  output  a  function  of  sun  position.  As  multiple  simulations 
for  different  sun  and  aircraft  positions  will  not  be  performed,  the  decision  to  not 
include  the  sun  is  chosen. 

The  overcast  day  provides  a  cloud  cover  which  absorbs  and  re-emits  radiation. 
A  cloudless  sky  should  be  a  selective  emitter  where  emission  bands  correspond  to 
the  atmospheric  constitutes,  for  example  water,  H20.  and  carbon  dioxide,  C 02 .  The 
denser  atmosphere  (cloud  cover)  fills  in  the  background  spectral  radiance  curve,  which 
better  approximates  a  blackbody,  as  it  has  more  particles  to  absorb  and  re-emit 
radiation.  Figure  3.27  shows  two  spectral  background  radiance  curves  produced  by 
MODTRAN,  of  a  typical  sunny  day  and  the  overcast  summer  day,  that  will  be  used 
in  the  simulation  for  the  MWIR  band.  Figure  3.28  displays  the  LWIR  case. 
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Figure  3.27:  MWIR  spectral  path  radiance  to  infinity  produced  by  MODTRAN  for 
(a)  a  sunny  day,  and  (b)  a  cloudy  day.  Note  how  the  presence  of  cloud  fills  in  the 
curve  to  approximate  a  greybody. 
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Figure  3.28:  LWIR  spectral  path  radiance  to  infinity  produced  by  MODTRAN  for 
(a)  a  sunny  day,  and  (b)  a  cloudy  day.  Note  how  the  presence  of  cloud  fills  in  the 
curve  to  approximate  a  greybody. 
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The  hemisphere  is  broken  into  ten  regions  in  elevation.  This  allows  different 
background  radiance  profiles  to  be  simulated  as  a  function  of  elevation.  Figure  3.29 
shows  the  integrated  background  radiance  as  a  function  of  elevation  from  the  earth 
surface  normal  where  90°  is  the  horizon.  Note  how  the  integrated  path  radiance  is 
a  function  of  elevation  angle  and  increases  as  the  path  moves  closer  to  the  horizon. 
This  is  because  the  path  becomes  more  optically  dense  as  more  atmosphere  at  low 
altitudes  is  included  in  the  simulation.  When  observing  at  zero  elevation  (zenith), 
the  optical  density  of  the  path  is  the  least,  thus,  the  path  radiance  curve  is  not  filled 
in  as  much  and  the  integrated  band  radiance  is  lower. 

3.5  Emission  and  Propagation  Algorithms 

During  a  simulation,  the  incident  and  reflected  radiance  of  a  facet  must  be  cal¬ 
culated.  A  tesselated  aircraft  model  allows  the  algorithm  to  be  facet-based  rather 
than  pixel-based  as  computer  graphics  algorithms  are.  Additionally,  a  pixel-based 
algorithm  (which  could  be  used  to  simulate  an  imaging  optical  system)  has  the  limi¬ 
tation  of  not  sampling  with  prior  knowledge  of  the  aircraft’s  geometry.  It  is  this  lack 
of  a  priori  information  that  increases  the  computation  time  of  blind  deterministic 
algorithms  and  drives  the  need  for  modified  algorithms  for  a  more  resource  friendly 
solution. 

A  detector  (whether  a  camera  or  missile  seeker-head)  that  is  orientated  in  the 
direction  of  the  aircraft  will  integrate  the  visible  facets  from  any  point  of  view  (POV) 
into  a  single  electrical  output  for  each  element  in  the  detector. 

As  discussed  in  Chapter  II,  the  emitted  radiation  from  a  surface  is  either  self- 
emitted  or  reflected.  The  algorithm  must  calculate  the  emitted  and  reflected  energy 
(in  radiance  or  intensity,  whichever  is  the  more  convenient  unit  for  further  use)  from 
each  facet.  To  express  this  mathematically  in  radiance  units,  Equation  (3.4)  lists 
the  observed  radiance  Te(observed)  as  a  sum  °f  the  self-emitted  radiance  and  reflected 
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Figure  3.29:  Integrated  path  radiance  to  infinity  produced  by  MODTRAN  for  (a) 
MWIR  and  (b)  LWIR  bands.  The  absolute  magnitude  of  the  values  is  not  important, 
simply  note  the  power  of  the  background  path  radiance  dependence  upon  elevation. 
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radiance. 
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Calculation  of  the  self  emitted  radiance  is  trivial.  The  radiant  thermal  emission  of 
Chapter  II  (at  some  surface  temperature,  T  [K ] )  is  scaled  by  the  specific  material 
emissivity  from  Equation  (2.4),  where  the  reflection  coefficient,  p(0),  is  taken  from 
the  material  DHR  measurements. 
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Conversely,  calculating  the  reflected  radiance  is  considerably  more  complex  and  the 
following  subsections  deal  with  various  methods  of  achieving  an  approximation  of 
the  analytical  solution.  The  accuracy  of  the  approximation  is  driven  chiefly  by  the 
sampling  resolution  of  the  model  and  functions  in  the  simulation. 

Analytically,  the  reflected  radiance  from  a  point,  x ,  in  the  direction  of  the 
observer,  is  the  hemispherical  integration  of  the  irradiance  at  a  point,  x,  and  the 
hemispherical  BRDF.  This  method  is  expressed  mathematically  in  Equation  (3.6), 
and  Figure  3.30  also  shows  this  graphically. 
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By  combining  Equations  (3.5)  and  (3.6),  an  accurate  definition  of  emitted  radiance 
results, 
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Observer  (Detector) 


facet  A 

Figure  3.30:  An  illustration  of  the  process  of  hemispherically 
integrating  the  irradiance  on  facet  A,  to  calculate  the  reflected 
radiance  in  the  direction  of  the  observer,  from  a  point,  x. 
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3.5.1  Specular  Reflection  Assumption.  The  simplest  way  of  evaluating  the 
reflected  radiance  integral,  Equation  (3.6),  is  that  of  a  purely  the  specular  surface, 
which  was  the  major  assumption  in  Bortle’s  and  Martinez’s  work,  [3]  and  [17].  This 
assumption  is  valid  when  a  highly  specular  surface  is  used;  however,  when  a  more 
diffuse  surface  is  used,  the  assumption  breaks  down  and  gives  inaccurate  results. 
Under  the  specular  assumption,  the  integral  is  not  required  to  be  evaluated  as  the 
BRDF  is  only  non-zero  when  the  incident  and  reflected  elevation  angles  are  equal  and 
the  azimuthal  angles  differ  by  180  degrees.  The  reflected  radiance  then  becomes  the 
incident  radiance  multiplied  by  the  reflectivity,  p{6). 
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This  method  is  very  fast  to  compute  as  only  a  single  incident  radiance  value  is  required. 


3.5.2  Radiosity.  The  first  algorithm  implemented  to  attempt  to  hemispher- 
ically  integrate  the  irradiance  contribution  was  classical  radiosity  as  described  in  [10] . 
Radiosity  is  an  algorithm  that  calculates  in  the  direction  of  flux  propagation  with  the 
intent  of  calculating  the  exitant  radiance  distribution  of  each  facet.  When  the  process 
has  finished,  a  value  of  observed  radiance  from  any  direction  can  simply  be  read  from 
the  data  file.  This  method  would  make  trend  analysis  very  easy,  simply  by  rotating 
around  the  model  in  azimuth  and  reading  off  radiance  values  from  each  facet. 

The  algorithm  starts  with  the  first  facet’s  self  emission,  where  the  isotropic 
radiance  distribution  is  calculated  from  Equation  (3.5),  as  a  function  of  p{0).  The 
facet’s  hemisphere  must  then  be  sampled  to  reveal  the  surrounding  geometry  (whether 
it  is  another  aircraft  facet  or  background).  This  method  divides  the  hemisphere  up 
into  a  2D  (0  and  6)  grid  where  the  exitant  radiance  is  assumed  to  be  constant  through 
the  sampled  solid  angle.  The  number  of  divisions  is  up  to  the  user,  however,  larger 
numbers  will  take  longer  to  calculate,  yet,  will  yield  a  more  accurate  solution.  This 
sampling  method,  more  commonly  known  as  ray  tracing,  is  described  in  detail  in 
Appendix  C. 
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For  solid  angles  where  background  is  visible  (a  possible  location  of  the  observer), 
the  emitted  radiation  is  stored  in  that  facet’s  observed  radiance  map  for  later  recall 
when  rendering  the  scene.  For  the  other  solid  angles  where  model  facets  are  visible, 
the  flux  to  each  facet  is  summed  and  the  algorithm  moves  to  the  next  step  for  each 
visible  facet.  Incident  flux  on  a  facet  can  be  easily  converted  to  irradiance  by  dividing 
by  the  facets  area.  From  there,  the  BRDF  is  scaled  by  this  irradiance  which  becomes 
the  exitant  radiance  from  that  facet.  The  process  of  sampling  the  hemispherical 
geometry  continues  as  before  and  the  algorithm  quickly  becomes  nested. 

The  algorithm  stops  when  the  self-emitted  radiance  from  the  first  facet  exits 
the  model  and  reaches  the  background  where  it  is  stored  in  the  observed  radiance 
map.  For  complex  models,  this  process  may  take  considerable  time;  thus,  stopping 
conditions  should  be  developed  to  truncate  the  random  walk  of  the  algorithm  through 
the  facet’s  geometry.  Once  the  algorithm  has  finished  propagating  the  first  facet’s  self 
emission,  it  moves  to  the  second  facet  and  restarts.  Stopping  conditions  for  the 
algorithm  can  include  truncating  the  random  walk  after  a  certain  number  of  bounces, 
or  ignoring  small  incident  flux  values  that  when  reflected  into  the  facet’s  hemisphere 
will  not  meaningfully  contribute  to  the  observed  radiance.  Figure  3.31  is  a  flow  chart 
that  outlines  the  process  in  graphical  form,  and  Equation  (3.9)  describes  the  process 
mathematically.  For  greater  detail  on  the  algorithm,  including  techniques  to  increase 
accuracy,  see  [10]  and  [5,  Chapter  6]. 
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Radiosity  algorithms  are  historically  slow,  as  path  lengths  carrying  significant  flux 
can  become  long.  Also,  they  are  suited  well  to  diffuse  reflections  where  the  sampling 
rate  can  be  kept  low.  However,  when  more  specular  BRDFs  are  used,  the  sampling 
issue  described  in  Section  3.3.2  becomes  important.  When  the  incident  irradiance  is 
multiplied  by  the  BRDF,  the  BRDF  must  be  adequately  described  (sampled).  This 
method  was  used  early  on  during  the  simulations  included  in  Appendix  A  and  provided 
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Figure  3.31:  The  steps  involved  in  the  Radiosity  algorithm  are  listed  above.  The 
number  of  bounces  becomes  important,  as  they  are  calculated  independent  of  the 
observation  direction.  Even  if  only  one  observation  is  required,  the  algorithm  must 
be  completed  in  full. 
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accurate  solutions.  However,  when  a  realistic,  specular  BRDF  was  introduced,  the 
sampling  rate  had  to  be  increased  dramatically.  This  increase  made  computation 
times  very  long,  even  for  the  cylinder  model.  For  this  reason,  radiosity  was  not  chosen 
as  an  algorithm  capable  of  simulating  IR  signatures  with  highly  specular  BRDFs. 

3.5.3  Point- Of- View  Hemispherical  Integration.  The  classical  radiosity  al¬ 
gorithm  of  [10]  has  an  easy-to-calculate-and-program  algorithm.  Because  of  this  prop¬ 
erty,  an  attempt  was  made  to  reverse  this  algorithm  starting  from  the  observer,  in 
the  hope  of  reducing  the  number  of  paths  and  computation  time. 

This  algorithm  starts  with  the  first  visible  facet  from  the  observer’s  POV.  From 
there,  Equation  (3.6)  is  approximated  by  evenly  sampling  the  hemisphere,  identically 
to  the  radiosity  case,  and  calculating  incident  radiance  on  the  facet.  To  include 
multiple  bounces,  by  calculating  the  reflected  radiance  of  Equation  (3.6)  again,  the 
algorithm  steps  into  the  visible  facets.  This  method  attempts  to  integrate  the  BRDF 
and  incident  radiance  using  the  method  described  in  Section  3.3  where  there  is  no 
a  priori  knowledge  of  the  reflectance  distribution.  The  same  issues  of  computation 
time  are  experienced,  and  this  blind  integration  method  is  rejected  as  a  suitable 
propagation  algorithm. 

3.5.4  Point- Of- View  Monte  Carlo  Hemispherical  Integration.  Similarly  to 
the  narrative  in  Section  3.3,  Monte  Carlo  hemispherical  integration  was  evaluated 
as  an  algorithm.  However,  as  the  BRDF  could  not  be  integrated  efficiently  with  a 
normal  or  uniform  distribution,  this  method  is  not  followed.  If  the  geometry  of  the 
simulation  created  many  paths,  for  example,  a  complex  cavity  like  an  engine  exhaust 
system,  then  this  method  may  prove  more  efficient  than  numerical  quadrature  as  the 
number  of  dimensions  and  bounces  becomes  large. 

3.5.5  Point- Of- View  Stratified  Hemispherical  Integration.  Section  3.3.3 
showed  that  the  most  efficient  method  of  integrating  the  BRDF,  one  of  the  terms 
in  Equation  (3.6),  is  to  use  stratified  numerical  quadrature  in  two  regions;  the  spec- 
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ular  peak  and  outer  diffuse  floor.  As  this  method  is  the  most  efficient  in  integrating 
the  BRDF,  it  foilows  that  it  will  be  the  most  efficient  in  calculating  Equation  (3.6). 
Hence,  the  best  propagation  algorithm  will  be  POV  stratified  hemispherical  integra¬ 
tion.  Figure  3.32  is  a  flow  chart  outlining  the  propagation  algorithm. 

This  algorithm  is  again  facet-based,  and  starts  at  the  first  facet  in  the  model. 
The  self-emitted  radiance  is  easily  calculated  using  Equation  (3.5)  and  the  reflected 
radiance  is  calculated  using  the  algorithm  in  Figure  3.32  to  approximate  the  integral 
in  Equation  (3.6).  Although  this  method  is  POV-based  and  has  to  be  recalculated 
for  each  observing  location,  an  answer  is  achieved  much  faster  than  in  the  classical 
radiosity  case. 

The  length  of  possible  paths  that  the  algorithm  can  calculate  when  integrating 
the  BRDF  over  the  whole  hemisphere  is  infinite.  The  infinite  length  paths  are  cre¬ 
ated  by  recurring  bounces  between  two  mutually  visible  facets.  The  diffuse  reflection 
calculations,  for  a  highly  specular  surface,  quickly  become  negligible,  even  after  the 
first  bounce. 

Two  approximations  can  be  made  to  reduce  the  number  of  paths.  The  first  is  to 
reduce  the  angular  search  space  of  each  bounce.  From  Equation  (3.6),  the  irradiance 
from  the  whole  hemisphere  must  be  integrated  to  produce  the  reflected  radiance.  How¬ 
ever,  if  the  surface  is  highly  specular,  then  the  angular  region  of  integration,  Q,  can 
be  reduced  and  still  calculate  an  accurate  result.  This  was  hinted  at  in  Section  3.3.3, 
however,  the  outer  region  integral  was  maintained.  Figures  3.19  and  3.20  show  the 
integral  values  over  the  peak  (0.917)  and  outer  (0.02)  regions,  respectively.  The  large 
value  for  the  peak  integral  compared  to  the  outer  result  proves  that  the  surface  (bare 
aluminum  in  this  case)  is  highly  specular.  Under  the  assumption  that  the  model  has 
a  highly  specular  surface  treatment,  then  the  solid  angle  can  be  reduced  to  only  the 
peak  region  of  the  BRDF,  with  a  negligible  effect  on  accuracy. 

The  second  approximation  is  to  stop  the  path  when  a  surface  of  very  high  emis- 
sivity  is  encountered.  In  this  case,  the  self-emitted  term  of  Equation  (3.4)  dominates, 
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Figure  3.32:  The  process  flow  for  the  hemispherical  numerical  quadrature  integra¬ 
tion  algorithm  is  less  complex  than  the  radiosity  one.  This  algorithm  also  allows  a 
solution  to  be  quickly  obtained  where  only  the  radiance  in  one  direction  is  calculated 
for  each  simulation. 
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rather  than  the  reflected  term;  when  this  is  the  case,  there  is  not  much  gained  from 
calculating  further  paths  as  the  result  will  be  scaled  by  the  very  small  reflectance. 
This  method  is  not  employed  in  the  algorithms,  however,  in  the  IR  range,  the  walls 
are  highly  emissive  (96%)  and  as  the  algorithm  does  not  calculate  a  bounce  off  the 
background  facets,  by  proxy,  this  method  was  implemented,  although,  it  was  not  for 
the  highly  emissive  model  facets. 

3.5.6  Memory  Considerations.  With  all  of  these  algorithms,  memory  allo¬ 
cation  cannot  be  ignored.  The  radiosity  algorithm  is  the  most  memory  intensive  as  it 
stores  the  emitted  radiance  maps  for  each  facet.  A  single  map  contains  a  value  for  each 
combination  of  9  and  </>.  The  cylinder  model  has  around  2000  facets,  and  typical  hemi¬ 
spherical  resolutions  were  200  in  each  dimension.  This  requires  2000  x  200  x  200  x  8 
bytes  (625Mb)  of  data.  To  have  this  much  memory  reserved  for  one  variable,  the 
computer  must  have  enough  memory  for  Mat  lab®  and  the  operating  system  to  still 
function,  without  using  swap  space  on  the  hard  drive.  If  swap  space  is  used,  then  all 
the  techniques  listed  here  to  save  processing  time  will  be  wasted. 

3.6  Artificial  Cavity  Blackbody  Source  Simulation 

When  calibrating  laboratory  systems,  a  perfect  blackbody  source  is  usually  re¬ 
quired.  However,  no  material  by  itself  exhibits  the  properties  of  a  perfect  blackbody 
as  there  is  always  a  reflected  component  as  the  emissivity  is  never  exactly  one.  Flat 
black  paint  is  very  close,  yet  it  is  still  not  accurate  enough  to  be  used  as  a  calibration 
source. 

A  useful  thermodynamics  property  is  that  radiation  inside  an  isothermal,  en¬ 
closed  cavity  matches  that  of  a  perfect  blackbody.  The  radiation  emitted  inside  the 
cavity  will  bounce  around,  being  absorbed  and  re-emitted,  until  the  radiation  inside 
the  cavity  is  spectrally  that  of  a  perfect  blackbody  source  at  that  same  temperature, 
even  when  very  low  emissivity  internal  surfaces  are  used.  To  use  this  radiation  as  a 
source,  a  small  hole  must  be  made  to  allow  the  radiation  to  escape.  This  hole  is  then 
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Table  3.3:  Variable  Definitions  for  the  Method  of 

Gouffe. 


Variable 

Definition 

£o 

Total  emissivity  of  cavity  source 

Emissivity  of  cavity 

£ 

Emissivity  of  internal  cavity  walls 

S 

Area  of  aperture 

s 

Surface  area  of  cavity  including  aperture 

So 

Surface  area  of  a  sphere  of  the 
same  depth  as  the  cavity 

l 

Depth  of  cavity 

r 

Radius  of  aperture 

emitting  perfect  blackbody  radiation.  The  radiation  becomes  less  representative  of  a 
perfect  blackbody  as  the  hole  becomes  larger.  There  are  many  theories  that  predict 
this  behavior,  however,  the  Method  of  Gouffe  is  provided  in  [11], 

Gouffe  predicts  the  total  emissivity,  £q,  from  a  cavity  source  as  a  function  of 
internal  emissivity,  e,  cavity  depth,  l,  hole  radius,  r,  area  of  cavity  aperture,  s,  internal 
surface  area  of  the  cavity,  S,  and  the  surface  area  of  a  perfect  sphere  with  the  same 
diameter  as  the  cavity,  So.  The  method  is  detailed  in  Equations  (3.10)  through  (3.12), 
and  the  variables  are  described  in  Table  3.3. 


So  —  £0  (1  +  k) 

(3.10} 

k  =  (1  —  e) 

s  s 

_s~s~0_ 

(3.11) 

/  S 

(3.12) 

°  ®  [1  -  5 

]  +  (f ) 

Equation  (3.12)  is  plotted  as  a  function  of  s/S  for  varying  internal  emissivities, 
c,  in  Figure  3.33,  which  is  reproduced  from  [26,  Figure  2-1].  The  ratio  of  l/r  is 
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Figure  3.33:  The  emissivity  of  the  blackbody  cavity  is  calculated,  using  Equa¬ 

tion  (3.12),  as  a  function  of  the  s/S  ratio  of  the  cavity  for  many  different  internal 
greybody  emissivities,  e.  The  lower  half  of  the  figure  calculates  the  s/S  ratio  as 
a  function  of  l/r  for  spherical,  conical  and  cylindrical  cavity  shapes.  Reproduced 
from  [26,  Figure  2-1]. 
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converted  to  the  s/S  ratio,  transferred  up  to  the  greybody  emissivity  of  the  cavity 
and  across  to  the  cavity  emissivity.  This  value  is  then  modified  to  account  for  cavity 
shape  by  k  in  Equation  (3.10).  The  different  aperture  sizes  create  different  s/S  ratios, 
and  thus,  eight  data  points  on  each  curve. 

The  curves  of  Figure  3.33  will  be  reproduced  by  the  model  in  Figure  4.1. 

3.7  A340-300  Simulation  Output 

The  output  of  the  simulation  is  a  spectral  intensity  curve  in  the  direction  of  the 
observer,  for  the  spectral  band  of  interest.  The  simulation  calculates  the  intensity 
at  the  target.  Thus,  to  calculate  the  apparent  intensity  at  the  detector,  a  spectral 
transmission  curve  must  be  applied.  To  gain  a  total  in-band  scalar  intensity  value, 
the  spectral  intensity  can  simply  be  integrated  across  the  waveband. 

To  calculate  the  intensity  in  the  direction  of  the  observer,  the  simulation  must 
convert  the  emitted  and  reflected  radiance  from  each  facet  to  intensity  by  summing 
the  radiances  and  multiplying  by  the  projected  facet  area.  This  process  is  expressed 
in  Equation  (3.13)  where  N  is  the  number  of  facets  in  the  model,  AjCos(^)  is  the 
projected  area  of  the  facet,  and  V)  is  the  visibility  function.  The  visibility  function 
is  a  binary  function:  1  if  the  facet  is  visible  from  the  view  point,  and  0  if  not.  The 
simulation  must  calculate  which  facets  are  visible  from  the  POV  and  apply  a  mask, 
V(i),  to  the  intensity  values  calculated  previously. 

N 

le  Target  ^  ^  (-^^  Self— Emitted  T  Tf'Refi  looted)  '  Aj  •  COS  (h; )  •  V) 

2=1 

3.8  Exercising  the  Model 

A  scenario  typical  of  an  aircraft  landing  at  an  international  airport  will  be 
simulated  to  calculate  the  apparent  irradiance  at  the  detector.  Then,  a  maximum 
value  for  noise  equivalent  irradiance  (NEI)  will  be  chosen  to  allow  a  system  to  track 
and  lock- on  to  the  aircraft.  A  lock-on  event  is  defined  as  when  the  apparent  irradiance 
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from  the  target  achieves  an  SNR  of  10,  or  ten  times  the  NEI  of  the  detector  system. 
An  expression  for  the  maximum  allowable  NEI,  and  yet  still  lock  onto  the  aircraft, 
is  developed  in  Equation  (3.14).  An  SNR  of  10  was  chosen  to  give  the  detector  a 
signature  large  enough  to  ensure  the  system  can  detect  and  track  the  aircraft. 
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The  NEI  is  a  parameter  that  describes  the  noise  power  of  the  entire  electro-optical 
system  of  the  detector.  This  allows  comparison  between  platforms  independent  of 
detector  type,  tracking  methods  and  noise  sources.  The  NEI  is  a  noise  power  param¬ 
eter  which  translates  the  internal  noise  sources  (e.g.,  Johnson,  shot  and  generation- 
recombination  noise)  as  an  equivalent  irradiance  at  the  optic. 


3. 8. 1  Contrast  Detection.  The  signal  detected  by  the  detector  is  the  sum  of 
the  target  and  background  irradiance  at  the  optic  plane  of  the  detector.  The  back¬ 
ground  exists  constantly  and  the  target  must  provide  contrast  against  the  background 
for  detection.  The  simulation  assumes  a  single-pixel  detector  where  the  whole  FOV 
is  integrated  into  one  signal. 

Figure  3.34  shows  the  irradiance  with  and  without  the  target  present.  The  first 
irradiance  signal  level,  Ee  NoTarget •  is  the  background  without  a  target  present. 


Ee  NoTarget 


=  E 


eBackground  DpQy 


(3.15) 


This  background  irradiance  is  calculated  with  MODTRAN  as  was  shown  in 
Figures  3.27  and  3.28.  The  integrated  background  radiance  across  the  band  then 
becomes  the  magnitude  of  the  signal.  The  ripples  on  the  signal  depict  the  electronic 
noise  that  is  present  in  all  systems,  such  as  Johnson  or  shot  noise.  The  NEI  represents 
the  power  contained  in  this  electronic  noise.  Thus,  for  the  target  to  be  detectable,  it 
must  have  a  signal  that  is  greater  than  the  random  noise.  The  signals  in  Figure  3.34 
are  expressed  as  irradiance.  Irradiance  is  calculated  by  Equation  (3.16)  where  the 
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Figure  3.34:  The  irradiance  signal  at  the  detector’s  optic  is  a 
combination  of  background  and  target  emissions. 

radiance  of  the  source  is  multiplied  by  the  projected  area  of  the  source  and  divided 
by  the  range  squared,  which  is  also  equal  to  the  radiance  multiplied  by  the  solid  angle 
that  the  source  occupies. 


rp  LeAs  cos(6»s)  to  _  1 e 

Ee  - - W - “  LSIfov  -  w 


(3.16) 


When  the  target  is  present  in  the  FOV,  the  background  irradiance  is  still  present, 
however,  the  target  is  now  blocking  some  solid  angle  of  background  and  must  be  sub¬ 
tracted  before  the  addition  of  the  signal  from  the  target.  To  calculate  the  propagated 
spectrum,  the  signal  from  the  target  must  be  multiplied  by  the  spectral  transmis¬ 
sion  discussed  in  Section  2.4.  Additionally,  the  path  radiance  from  the  target  to  the 
detector  must  be  added  to  the  signal.  This  path  radiance  replaces  the  subtracted 
background  radiance  over  the  solid  angle  of  the  target.  The  solid  angle  of  the  target 
is  calculated  by, 

^Target  =  (3-17) 
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Combining  these  terms  gives  the  irradiance  with  the  target  present 


Ee 
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(3.18) 


As  introduced  above,  the  detection  process  operates  on  contrast.  A  threshold  is 
set  above  the  background  irradiance  level  and  when  a  target  is  present  in  the  FOV,  the 
signal  from  the  target  must  lift  the  observed  signal  over  the  threshold.  A  threshold  of 
ten  times  the  noise  level  is  chosen  to  guarantee  not  only  detection  but  enough  signal 
power  for  tracking  algorithms  to  function  correctly.  This  SNR  of  10  is  arbitrary,  where 
the  actual  required  value  depends  upon  the  sophistication  of  detection  and  tracking. 

Because  the  detection  process  is  only  concerned  with  the  difference,  the  back¬ 
ground  signal  can  be  normalized  out  and  the  detection  signal  as  shown  in  Figure  3.34 
is  the  difference  of  the  irradiance  with  and  without  the  target, 
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^Background  ^Target 


(3.19) 


3.8.2  Flight  Path.  A  standard  Instrument  Landing  System  (ILS)  approach 
path  is  simulated.  ILS  approaches  are  directed  by  a  radio  frequency  (RF)  system 
where  the  aircraft  navigates  down  a  virtual  3°  glideslope,  starting  about  four  nautical 
miles  (7km)  out  from  the  runway.  The  aircraft  will  also  be  given  an  angle  of  attack 
of  +10°,  which  is  approximately  representative  of  the  attitude  of  an  aircraft  on  final 
approach. 

3. 8. 3  Aircraft  Engagement  Positions.  Three  positions  are  chosen  to  simulate 
the  engagement  of  a  MANPADS  system  against  the  aircraft.  The  first  position  is  1km 
from  the  threshold  at  an  altitude  of  50m,  the  second  is  at  2km  out  at  100m  altitude, 
and  the  third  is  3km  out  and  150m  in  altitude.  The  detector/observer  is  located  2km 
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Figure  3.35:  Three  engagement  positions  are  chosen  to  calcu¬ 
late  the  infrared  signature  irradiance  at  the  detector. 

from  the  threshold  and  500m  off  the  centerline  of  the  runway.  These  positions  are 
illustrated  in  Figure  3.35. 

The  positions  were  chosen  as  examples  only  and  are  not  meant  to  represent  a 
complete  parametric  study  of  each  of  the  variables.  They  are  merely  a  demonstration 
of  the  application  of  an  IR  signature  simulation  code. 

3.9  Chapter  Summary 

This  Chapter  has  applied  the  theory  introduced  in  Chapter  II  to  the  methods, 
models  and  algorithms  that  will  be  used  in  the  IR  signature  simulations.  The  results 
of  these  simulations  are  detailed  and  discussed  in  Chapter  IV. 
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IV.  Results  and  Analysis 


This  chapter  details  the  results  of  the  different  simulations  performed  using  the 
algorithms  developed  in  Chapter  III.  The  aim  of  this  thesis  is  to  demonstrate 
the  effect  of  the  BRDF  and  multiple  bounce  calculations  in  IR  signature  simulations. 
Several  different  simulations  are  performed  to  aid  in  the  demonstration  of  the  ob¬ 
jectives.  Initially,  simple  geometric  cases  (where  analytical  results  can  be  obtained) 
were  simulated  to  ensure  the  software  was  coded  correctly,  and  that  the  propagation 
algorithm  is  a  valid  approximation  to  the  reflected  radiance  integral  of  Equation  (3.6). 
The  results  of  these  simple  simulations  can  be  found  in  Appendix  A,  where  the  radios- 
ity  algorithm  performed  quite  well  on  diffuse  surfaces.  However,  it  failed  to  calculate 
the  correct  results  for  specular  surfaces  in  a  timely  manner. 

4-1  Artificial  Cavity  Blackbody  Source  Simulations 

A  more  complex  simulation  than  the  simple  analytical  proofs  was  performed  to 
check  the  multiple  bounce  algorithm,  and  also  as  a  quantitative  check  of  the  radio- 
metric  calculations  in  the  algorithm. 

This  property  is  simulated  to  test  the  accuracy  of  the  algorithm,  using  several 
facetised  spherical  cavities  with  varying  aperture  sizes.  Calculations  of  s/S  ratios  for 
the  other  aperture  sizes,  and  a  sample  simulated  emissivity  calculation  are  included 
in  Appendix  B.  Figure  4.1  attempts  to  simulate  Figure  3.33  for  the  e  =  0.5  case. 
Note  the  effect  of  the  number  of  bounces  in  the  simulation.  After  only  one  bounce, 
the  cavity  emissivity  is  approximately  75%  of  the  theoretical  result.  The  simulated 
emissivity  of  the  cavity  increases  after  each  further  bounce.  After  four  bounces,  the 
simulation  is  96%  accurate;  the  ability  of  the  propagation  algorithm  to  accurately 
predict  the  analytical  emissivity  is  shown.  Further  bounces  were  not  calculated  due 
to  processing  time,  and  also  because  the  algorithm  was  predicting  the  cavity  emissivity 
correctly.  As  the  reflectance  distribution  of  the  cavity  is  diffuse,  the  whole  hemisphere 
contributes  equally  with  respect  to  the  BRDF,  thus,  the  whole  hemisphere  of  each 
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Figure  4.1:  The  emissivity  of  the  e  —  0.5  greybody  cavity 

is  simulated  using  the  propagation  algorithm  of  Equation  (3.7). 

The  upper  curve  is  the  analytical  limit  (Equation  3.12)  of  the 
cavities’  ability  to  approximate  a  perfect  blackbody  as  a  func¬ 
tion  of  s/S.  Several  different  aperture  sizes  were  simulated  to 
produce  a  curve  in  stead  of  a  single  point.  After  each  bounce 
the  apparent  cavity  emissivity  increases  as  the  approximation  of 
the  reflected  radiance  integral  is  more  precise. 

facet  must  be  sampled.  This  creates  a  large  number  of  paths  of  appreciable  radiance 
that  must  be  calculated  to  obtain  the  correct  result. 

In  reality,  the  greybody  radiation  in  the  cavity  has  experienced  an  almost  lim¬ 
itless  number  of  bounces  as  radiation  is  reflected,  absorbed  and  emitted.  Referring 
the  simulation  back  to  Equation  (3.4),  the  self-emitted  radiance  of  the  visible  facets 
inside  the  cavity  is  attenuated  by  the  greybody  emissivity,  z  =  0.5.  The  reflected 
component,  then  becomes  considerable  as  the  number  of  bounces  increases,  the  paths 
become  longer  and  the  irradiance  on  the  visible  facets  increases  to  provide  the  other 
half  of  the  blackbody  spectrum. 
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4-2  IR  Range  Simulations 

Bortle  [3]  performed  experiments  in  AFRL’s  IR  range  to  validate  the  LCAIR 
code  developed  by  Martinez  [17].  Detailed  work  was  performed  to  produce  calibrated 
and  normalized  radiance  images  of  the  cylinder  model  in  the  range.  Images  were  taken 
in  the  MWIR  and  LWIR  spectral  regions  and  compared  to  the  output  of  LCAIR.  In 
each  of  the  figures  presented  for  the  IR  range  simulations,  the  physical  range  image 
will  be  presented  beside  the  simulated  image  for  comparison. 

The  following  sections  will  detail  limitations  and  advantages  of  the  single  spec¬ 
ular  bounce  LCAIR  algorithm,  the  multi  bounce  POV  stratified  hemispherical  inte¬ 
gration  algorithm,  and  combinations  of  both. 

4-2.1  Single  Bounce  and  Specular  Reflectance.  To  begin,  Bortle  found  that 
the  single  bounce,  specular  reflectance  LCAIR  algorithm  performs  well  when  the 
model’s  surface  is  highly  specular  and  the  specular  bounce  hits  a  highly  emissive 
facet,  whether  model  or  background. 

The  gloss  black  painted  model  is  simulated  the  best  by  LCAIR.  Although  it  is 
highly  specular,  as  can  be  seen  in  the  peak  of  the  BRDF  (Figure  2.12),  it  is  also  rea¬ 
sonably  emissive  over  all  wavelengths  which  was  seen  in  the  DHR  data  in  Figure  2.10. 
Figure  4.2  shows  the  simulation  of  the  gloss  black  painted  cylinder  model  in  the  IR 
range.  Visibly,  the  single  bounce  specular  reflectance  simulation  is  a  good  match  to 
the  physical  image.  This  is  achieved  because  the  surface  is  specular,  and  where  the 
specular  bounce  hits  other  model  facets,  the  emissivity  of  the  gloss  black  paint  is 
high  (see  Figure  2.10).  The  considerable  emissivity  of  the  gloss  black  paint  makes  the 
single  bounce  algorithm  approximate  the  Les; (Emitted)  term  of  Equation  (3.6)  much 
more  accurately.  Remember  that  the  observed  radiance  is  a  sum  of  the  emitted 
and  reflected  radiance,  and  the  emissivity  and  reflectance  coefficients  sum  to  one, 
e(9)  +  p{9 )  =  1.  Thus,  when  one  of  these  dominates,  the  other  term  can  be  considered 
negligible.  Conversely,  when  the  reflectance  term  dominates,  the  self-emitted  com- 
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(a)  Single  Specular  Bounce  (b)  IR  Range  Image 


Figure  4.2:  The  single  specular  bounce  assumption  works  well 
for  the  gloss  black  case  where  the  model  is  reasonably  emissive, 
where  all  specularly  reflected  rays  hit  an  emissive  facet,  and  the 
reflected  radiance  term  is  approximated  well.  Colormap  units 
are  [ W/(m 2  —  sr)].  Subfigure  (b)  is  reproduced  from  [3]. 

ponent  becomes  negligible  and  the  reflected  radiance  term  dominates  the  signature, 
similar  to  a  mirror. 

Although  the  emissivity  of  the  gloss  black  is  only  60-80%,  it  numerically  con¬ 
tributes  around  92%  of  the  final  result.  This  is  because  of  the  temperature  difference 
between  the  model  and  the  IR  range  walls.  The  model,  on  average,  is  300K,  and 
the  walls  260K.  Figure  4.3  plots  the  blackbody  spectral  radiance  curves  for  those  two 
temperatures.  The  vertical  lines  define  the  MWIR  (3  —  5/im)  and  LWIR  (8  —  12 /im) 
spectral  regions.  Note  that  the  300K  curve  is  much  higher  than  the  260K  curve,  and 
when  integrated  over  the  MWIR  band,  the  exitant  radiance  for  the  300K  blackbody 
is  1.92  \W/(m 2  —  sr)],  and  0.37  \W/{m 2  —  sr)]  for  the  260K  blackbody.  If  the  emis¬ 
sivity  and  reflectance  were  equal,  then  the  model  self  emission  (T  =  300K)  is  still 
significantly  larger  than  the  reflected  radiance  from  the  background.  However,  in  the 
gloss  black  case,  the  emissivity  is  60-80%  and  this  effect  is  magnified,  where  the  self 
emitted  radiance  is  1.3  [W/(m2  —  sr)\  and  the  reflected  radiance  is  0.11  [W/  (m2  —  sr)]. 
This  is  the  reason  that  the  algorithm  can  use  a  single  bounce  algorithm  and  still  pro- 
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Figure  4.3:  The  spectral  radiance  of  a  260K  and  300K  black- 
body  source  are  shown.  The  vertical  lines  show  the  extents  or 
the  two  spectral  bands,  MWIR  (3  —  5 pm)  and  LWIR  (8  —  12 fim). 

Colormap  units  are  [W/(m2  —  sr)]. 

duce  suitable  results  (over  90%  correct).  The  reflected  ray  will  always  hit  an  emissive 
surface,  whether  background  or  another  model  facet. 

However,  the  single  bounce  assumption  fails  in  the  bare  aluminum  case,  where 
the  model  is  highly  (approximately  90%)  reflective.  In  this  case,  the  reflected  term 
dominates  the  observed  radiance.  Even  with  the  MWIR  integrated  radiance  from 
Figure  4.3,  the  reflected  radiance  is  now  0.33  \W/{m2  —  sr)],  and  the  self-emitted 
radiance  is  reduced  to  0.19  [W/(m2  —  sr)].  Thus,  with  the  single  bounce,  where  this 
reflected  radiance  is  dropped,  the  result  underpredicts  by  60%.  Figure  4.4  shows 
this  effect;  the  black  regions  are  where  the  vector  from  the  observer  to  the  facet  is 
specularly  reflected  onto  one  of  the  wing  facets.  This  has  caused  the  observed  radiance 
from  these  facets  to  be  severely  under  predicted,  through  the  process  described  above. 
The  model  is  rotated  to  a  pitch  of  +45°  and  this  is  why  the  top  of  the  cylinder  is 
visible.  It  is  also  why  the  black  region  is  not  only  in  line  with  the  wings  in  the  vertical 
dimension. 
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(c)  Single  Specular  Bounce  (70°  azimuth)  (d)  IR  Range  Image  (70°  azimuth) 

Figure  4.4:  The  single  bounce  specular  assumption  fails  when  the  specularly  re¬ 

flected  ray  hits  a  highly  reflective  (cylinder  model  wing)  facet.  In  this  case,  the  re¬ 
flected  radiance  term  is  under  predicted  by  60%.  Colormap  units  are  [' W/(m 2  —  sr)]. 
Subfigures  (b)  and  (d)  are  reproduced  from  [3]. 
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The  error  is  caused  by  the  approximation  of  Equation  (3.6).  The  Le^  Emitted 
term  is  the  incident  radiance  on  facet  A  from  facet  B  (refer  to  Figure  3.30  for  facet 
definitions).  A  single  bounce  calculation  only  calculates  the  self-emitted  portion  of 
Equation  (3.4)  and  ignores  the  reflected  term.  The  error  occurs  when  the  reflectance 
of  facet  B  is  high.  Thus,  a  large  portion  of  emitted  radiance  is  lost  and  the  blackbody 
radiance  is  reduced  by  the  low  emissivity  value. 

A  second  bounce  is  required  to  evaluate  the  irradiance  on  facet  B  to  calculate 
the  reflected  radiance  integral  of  Equation  (3.6)  to  complete  the  Le^  Emitted  term  of 
facet  B.  This  process  is  nested,  and  many  bounces  may  have  to  be  calculated  (through 
multiple  bounces)  to  properly  evaluate  the  Le^  Emitted  term  for  reflective  surfaces.  In 
theory,  an  infinite  number  of  bounces  are  required,  however,  several  assumptions  can 
be  made  to  truncate  the  iteration  and  provide  a  timely  solution. 

4-2.2  Multiple  Bounces  and  Specular  Reflectance.  To  improve  on  the  sim¬ 
plistic  LCAIR  algorithm  (a  single  specular  bounce) ,  multiple  bounces  were  included  in 
the  simulation,  whilst  still  keeping  the  specular  assumption.  Multiple  bounce  calcula¬ 
tions  do  not  drop  the  reflected  radiance  term  from  the  observed  radiance  equation,  as 
Equation  (3.6)  is  nested  inside  itself  for  each  specular  bounce.  Figure  4.5  shows  the 
results  when  two  specular  bounces  are  calculated  from  the  same  simulation  as  before 
(LWIR  bare  aluminum). 

The  second  bounce  has  added  to  the  emitted  radiance  of  the  second  facet.  In 
these  simple  cases,  the  second  bounce  hits  the  background  of  the  IR  range,  and  the 
algorithm  stops.  The  highly  reflective  aluminum  creates  the  need  to  keep  tracking 
bounces  throughout  the  model.  In  this  case,  only  two  bounces  are  required  for  all 
paths  to  exit  to  the  background,  however,  with  more  complicated  geometry,  the  num¬ 
ber  of  bounces  required  would  be  more,  and  will  be  explored  with  the  A340-300  model 
simulations. 

Another  demonstration  of  how  the  second  bounce  contributes  to  the  accuracy 
of  the  bare  aluminum  in  the  MWIR  band  simulation  is  shown  in  Figure  4.6.  In 
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(a)  Single  Specular  Bounce  (0°  azimuth)  (b)  Single  Specular  Bounce  (70°  azimuth) 


(c)  Two  Specular  Bounces  (0°  azimuth)  (d)  Two  Specular  Bounces  (70°  azimuth) 


(e)  IR  Range  Image  (0°  azimuth)  (f)  IR  Range  Image  (70°  azimuth) 


Figure  4.5:  By  including  a  second  specular  bounce,  all  reflected  rays  now  exit  to  a 
highly  emissive  background  facet.  Thus,  the  reflected  radiance  term  of  Equation  (3.6) 
is  approximated  well.  Colormap  units  are  [ W/(m 2  —  sr )].  Subfigures  (e)  and  (f)  are 
reproduced  from  [3]. 
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(e)  IR  Range  Image  (f)  IR  Range  Image 


Figure  4.6:  By  including  a  second  specular  bounce,  all  reflected  rays  now  exit  to  a 
highly  emissive  background  facet.  Thus,  the  reflected  radiance  term  of  Equation  (3.6) 
is  approximated  well.  Colormap  units  are  [W/(m2  —  sr)].  Subfigures  (e)  and  (f)  are 
reproduced  from  [3]. 
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Figure  4.6  (a),  where  the  first  bounce  hits  another  model  facet,  the  300K  self  emission 
is  offset  by  the  low  emissivity  of  the  aluminum  and  the  result  approximately  equals 
the  96%  emission  of  the  range  walls.  This  is  the  reason  why  there  are  no  vertical 
black  stripes  in  the  simulation,  and  the  models  appear  to  be  evenly  colored. 

This  is  purely  a  coincidence,  and  is  a  result  of  wave  band,  facet  temperature, 
reflectivity  and  colormap  choice.  When  the  second  bounce  is  added,  the  observed 
radiance  (of  facets  which  have  another  facet  as  the  specular  reflection)  is  increased 
beyond  the  others  and  produces  the  yellow  glow  that  can  be  seen  in  the  range  image, 
Figure  4.6  (e  and  f).  The  result  is  an  algorithm  that  is  a  closer  simulation  of  the 
physical  processes. 

For  purely  specular  reflections,  the  cylinder  model  only  requires  two  bounces 
for  all  sample  rays  to  exit  to  the  background.  After  the  two  bounces,  all  rays  meet 
the  highly  emissive  background  and  Equation  (3.6)  is  suitably  approximated. 

4-2.3  Single  Bounce  and  Integrated  BRDF.  Figures  4.5  (a  and  b)  show  the 
single  bounce  simulations  for  the  bare  aluminum  LWIR  simulations.  The  change  in 
magnitude  of  observed  radiance  between  single  and  double  bounce  paths  is  sharp. 
This  gradient  is  a  result  of  the  specular  assumption,  where  one  ray  is  used  to  sample 
(represent)  the  whole  hemisphere.  In  reality,  there  is  a  gradual  change  from  a  facet 
which  has  all  background  in  its  hemispherical  view  and  one  which  has  some  model 
facets.  Additionally,  the  specular  lobe  has  some  width  to  it  (recall  Figure  2.9),  thus, 
significant  reflected  radiance  can  come  from  off  specular  angles. 

Relaxing  the  strict  specular  assumption  and  integrating  across  the  BRDF  smooths 
this  gradient,  to  improve  the  algorithm  further.  Figure  4.7  shows  the  result  of  inte¬ 
grating  across  the  hemisphere  with  the  BRDF  inside  Equation  (3.6),  for  the  bare 
aluminum  LWIR  simulation,  with  a  single  bounce.  The  result  is  that  the  gradient 
between  radiance  values  has  been  smoothed,  but  the  single  bounce  reflected  radiance 
contribution  failure  is  still  present.  This  figure  was  included  to  show  that  just  purely 
integrating  over  the  BRDF  with  a  single  bounce  calculation  does  not  automatically 
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(a)  Single  Specular  Bounce 


(b)  Single  Integrated  BRDF  Bounce 


Figure  4.7:  (a)  Zoomed-in  image  of  Figure  4.5  (a).  Even  when  the  BRDF  is  inte¬ 

grated  over  a  single  bounce  (b),  the  same  issues  from  a  highly  reflective  facet  occur. 
This  indicates  that  the  specular  assumption  is  quite  valid  as  the  two  results  look  sim¬ 
ilar.  The  effect  of  the  integration  can  be  seen  in  the  gradient  between  the  yellow  and 
brown/black  facets.  In  (b),  there  is  a  smoother  gradient  as  there  are  some  red  facets 
that  separate  the  yellow  and  brown/black  ones.  Colormap  units  are  \W/(m2  —  sr)]. 
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provide  an  accurate  simulation;  integrating  the  BRDF  over  multiple  bounces  is  re¬ 
quired  for  the  highest  fidelity  solution. 

4.2.4  Multiple  Bounces  and  Integrated  BRDF.  To  gain  the  best  result  pos¬ 
sible,  multiple  bounces  whilst  integrating  the  BRDF  must  be  calculated.  Previously, 
the  cylinder  model  can  only  support  two  specular  bounces  before  the  paths  exit  to 
background.  However,  by  integrating  the  BRDF  over  the  full  hemisphere,  an  infinite 
number  of  paths  and  bounces  are  created.  Section  3.5.5  justified  the  choice  for  trun¬ 
cating  the  number  of  paths  (as  highly  specular  surfaces  are  being  simulated)  through 
only  integrating  across  the  BRDF  peak  region. 

Calculating  the  multiple  bounces  (whilst  integrating  the  BRDF)  takes  more 
time  than  using  the  purely  specular  assumption;  however,  as  shown  in  the  previous 
sections,  it  is  more  accurate.  Similar  to  the  multiple  specular  bounce  algorithm, 
Figure  4.8  shows  the  LWIR  bare  aluminum  result  when  two  bounces  are  integrated 
across  the  BRDF  peak  region,  previous  figures  are  included  for  reference  as  subfigures. 
The  bare  aluminum  MWIR  result  is  included  in  Figure  4.9  also  for  comparison  to 
Figure  4.6  where  two  specular  bounces  are  calculated. 

As  the  number  of  bounces  is  increased,  the  number  of  paths  also  increases, 
approximating  the  reflected  radiance  integral  more  accurately.  Figure  4.10  shows  the 
reflected  radiance  term  of  facet  ^1388  converging  to  the  final  result  as  the  number 
of  bounces  increases.  Previously,  only  two  specular  bounces  were  possible,  however, 
with  the  wider  angular  search  region,  up  to  four  bounces  are  possible.  These  further 
bounces  have  an  appreciable  effect  on  the  integral  approximation  as  they  are  in  the 
peak  region.  The  convergence  rate  slows  as  the  geometry  becomes  more  complex,  thus 
more  paths  (more  bounces)  are  required  to  evaluate  the  reflected  radiance  integral 
properly. 

4.2.5  Algorithm  of  Choice.  From  the  results  above  and  Chapter  III,  the 
most  efficient  algorithm  is  very  much  dependent  upon  the  BRDF  of  the  material.  The 
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(a)  Two  Specular  Bounces 


(b)  Two  BRDF  Bounces 
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Figure  4.8:  The  addition  of  the  second  integrated  BRDF  bounce  fills  in  the  under 
prediction,  similar  to  the  specular  reflectance  case.  When  compared  to  the  real  image, 
the  two  BRDF  bounce  algorithm  performs  well,  however,  it  is  not  superior  to  the 
double  bounce  specular  algorithm.  Colormap  units  are  [W/(m 2  —  sr)].  Subfigure  (c) 
is  reproduced  from  [3]. 


(c)  IR  Range  Image 
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(c)  IR  Range  Image 


Figure  4.9:  The  MWIR  result  is  similar  to  the  LWIR  simulation,  the  second  BRDF 
bounce  provides  the  increased  reflection  in  the  wing  fillet  area.  The  MWIR  range 
image  is  a  much  higher  resolution  image  than  the  LWIR  image.  This  increased  res¬ 
olution  allows  a  greater  comparison  of  the  simulation  results.  Colormap  units  are 
\W/{m2  —  sr)].  Subfigure  (c)  is  reproduced  from  [3]. 


97 


14 


Figure  4.10:  The  convergence  to  the  final  result  (to  match 

the  IR  Range  image)  for  facet  #1388  is  shown.  Note  that  it 
takes  four  bounces  to  achieve  the  correct  result  as  a  result  of 
integrating  the  BRDF  over  the  peak  region. 

computation  time  is  also  dependent  upon  the  number  of  paths  that  carry  appreciable 
flux;  for  diffuse  surfaces,  this  is  number  can  be  considerable,  and  for  more  specular 
surfaces,  the  number  is  reduced. 

The  A340-300  to  be  simulated  in  the  next  section  is  assumed  to  be  constructed 
with  bare  aluminum.  Thus,  the  specular  assumption  will  be  valid  and  is  used. 


4-3  A340-300  Simulations 

The  final  simulation  that  was  performed  is  the  engagement  scenario  of  the  A340- 
300  by  a  threat  located  outside  a  major  airport.  Section  3.7  dealt  with  the  setup 
and  theory  of  the  simulation,  and  the  results  are  included  below.  The  simulations  are 
performed  to  show  the  end  application  of  the  IR  simulation  process.  The  scenarios  are 
broken  into  MWIR  and  LWIR  simulations,  and  for  position  1,  the  effect  of  increasing 
number  of  bounces  will  also  be  shown.  The  process  of  calculating  the  result  will  be 
shown  for  position  1,  whilst  only  the  final  results  for  positions  2  and  3  will  be  shown. 
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Figure  4.11:  The  output  of  the  MWIR  simulation  at  Position 
1  is  target  intensity  with  no  atmospherics  included.  Note  the 
wavelength  of  the  peak  in  the  curve.  4 /im  corresponds  to  the 
peak  emission  of  a  724K  blackbody  source.  This  temperature  is 
a  weighted  combination  of  the  exhaust  plume  and  turbine  hot 
parts  annulus.  The  background  spectral  peak  is  dwarfed  by  the 
magnitude  of  the  hotter  Planckian  curve,  and  is  not  visible. 

4-3.1  Position  1  -  MWIR.  The  first  position  simulates  when  the  observer 
(detector)  is  behind  the  aircraft  (refer  to  Figure  3.35),  where  the  detector  is  1km 
behind,  500m  to  the  side  and  50m  below  the  aircraft.  The  aircraft  also  has  a  pitch  of 
10°  to  simulate  the  landing  attitude. 

The  MWIR  simulation  (3  —  5 /im)  is  performed  and  the  results  are  included 
below.  The  output  of  the  propagation  algorithm  is  spectral  intensity,  Ie  Target  from 
Equation  (3.13),  observed  from  the  wireframe  model  with  no  atmospherics  added. 
This  spectral  intensity  is  shown  in  Figure  4.11.  The  path  transmission,  as  computed 
by  MODTRAN,  is  applied  to  gain  the  apparent  target  intensity,  and  is  then  converted 
to  apparent  target  irradiance  at  the  detector  with  Equation  (3.16).  The  results  of 
this  process  are  shown  in  Figure  4.12.  Once  the  spectral  apparent  target  irradiance 
is  calculated,  the  background  and  path  radiances,  again  computed  by  MODTRAN, 
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Apparent  Target  Irradiance 


Figure  4.12:  Spectral  transmission,  as  computed  by  MODTRAN,  is  applied  to  the 
target  intensity,  from  Position  1,  and  is  converted  to  irradiance  at  the  optic. 


100 


3 


A  \fim\ 

Figure  4.13:  The  spectral  irradiance  as  measured  at  the  de¬ 

tector,  from  Position  1,  with  all  atmospherics  included. 

are  passed  to  Equation  (3.18)  to  gain  the  total  received  irradiance  at  the  optic  which 
is  shown  in  Figure  4.13. 

These  figures  all  show  spectral  irradiance  or  intensity.  To  calculate  the  output  of 
the  single-pixel  detector,  the  spectrum  must  be  integrated  to  calculate  total  irradiance 
across  the  band.  To  give  a  feel  for  the  scalar  magnitudes  involved  between  the  different 
terms  of  Equation  (3.18),  they  are  listed  in  Table  4.1  for  each  simulation.  The  total 
MWIR  irradiance  received  at  the  detector  is  3.034  x  10_3lF/m2;  thus,  the  required  NEI 
with  a  SNR  of  10  is  3.034  x  10_4lF/m2.  The  integrated  irradiance  in  the  MWIR  band 
does  not  vary  within  three  significant  figures  when  varying  the  number  of  calculated 
bounces.  This  is  because  the  observation  angle  only  supports  specularly  reflected 
paths  that  hit  the  exhaust  plume  or  hot  parts  annulus  on  the  first  bounce.  There 
are  paths  with  up  to  four  bounces  possible,  however,  they  do  not  carry  significant 
radiance  and  do  not  contribute  appreciably  to  the  result. 

Figure  4.14  shows  the  image  of  the  aircraft  rotated  to  the  view  angle  of  position  1. 
An  image  where  the  background  is  also  drawn  is  included  to  show  contrast  between 
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Table  4.1:  Position  1  Detector  Irradiance. 


Band 

# 

Bounces 

Ep  Tgt  App 

[W/m2] 

Ee  Bg  App 

[W/m2] 

Ep  Path  App 

[W/m2] 

Ep  Total  App 

[W/m2] 

MWIR 

1 

3.03  x  10"3 

5.64  x  10"7 

3.23  x  10"7 

3.03  x  10“3 

2 

3.03  x  10"3 

5.64  x  10"7 

3.23  x  10"7 

3.03  x  10“3 

3 

3.03  x  10“3 

5.64  x  10“7 

3.23  x  10“7 

3.03  x  10“3 

10 

3.03  x  10"3 

5.64  x  10"7 

3.23  x  10"7 

3.03  x  10“3 

LWIR 

10 

3.61  x  10“3 

6.55  x  10“5 

2.23  x  10“5 

3.56  x  10“3 

the  aircraft  and  background  irradiance.  It  is  important  to  remember  that  the  detector 
has  to  integrate  over  the  whole  FOV,  which  includes  the  background.  Figure  4.14  (b) 
shows  this  view,  where  the  exhaust  direct  emission  and  reflections  off  the  fuselage  are 
what  lifts  the  signature  out  of  the  atmospheric  background  and  electronic  noise  to 
enable  detection. 

4-3.2  Position  1  -  LWIR.  The  previous  section  demonstrated  the  process 
and  results  for  the  MWIR  case.  For  the  LWIR  case,  the  process  is  the  same,  however, 
the  balance  between  the  signature  elements  differs  and  magnitudes  change.  In  the 
LWIR  band,  the  background  becomes  more  important  as  the  background  self  emission 
peaks  in  the  LWIR  band.  Figure  4.15  shows  the  spectral  target  intensity  for  the  LWIR 
case,  and  is  applied  to  the  spectral  transmission  of  Figure  4.16  to  calculate  the  total 
apparent  irradiance  at  the  optic  from  Position  1,  which  is  shown  in  Figure  4.17. 
Figure  4.20  shows  the  view  of  the  aircraft  at  Position  1  in  this  LWIR  case.  Note  the 
increased  contribution  from  the  terrain.  The  terrain’s  temperature  is  300/7  which  has 
a  self-emission  peak  at  10 \xm.  This  peak  lies  in  the  middle  of  the  LWIR  band,  hence, 
the  increased  earth  shine  reflection  off  the  aircraft  compared  to  the  MWIR  case. 
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(a)  Aircraft  Only 
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(b)  Aircraft  plus  Background 

Figure  4.14:  The  two  MWIR  images  show  the  orientation  of  the  aircraft  at  Position 
1  when  viewed  from  the  observer  with  (a)  the  background  removed,  and  (b)  the 
background  included  to  show  contrast  between  the  aircraft  emission  and  reflection, 
and  the  background  irradiance.  Colormap  units  are  [ W/(m 2  —  sr)]. 
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Figure  4.15:  The  output  of  the  LWIR  simulation  from  Position  1  is  target  intensity 
with  no  atmospherics  included. 
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Path  Transmission 


Figure  4.16:  Spectral  transmission,  as  computed  by  MODTRAN,  is  applied  to  the 
target  intensity  and  converted  to  irradiance  at  the  optic  from  Position  1. 
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Figure  4.17:  The  total  irradiance  as  measured  at  the  optic,  from  Position  1,  with 
all  atmospherics  included. 
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(a)  Aircraft  Only 
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Figure  4.18:  The  two  LWIR  images  show  (a)  the  background  removed,  and  (b)  the 
background  included  to  show  contrast  between  the  aircraft  emission  and  reflection, 
and  the  background  irradiance  at  Position  1.  Note  the  increased  contribution  in  the 
LWIR  band  from  the  ground  reflection  off  the  aircraft.  The  ground  (302. 5K)  has  its 
peak  emission  at  about  10 /im.  Colormap  units  are  [ W/(m 2  —  sr)]. 


(b)  Aircraft  plus  Background 
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Table  4.2:  Position  2  Detector  Irradiance. 


Band 

# 

Ef.  Tgt  App 

Ee  Bg  App 

Ef.  Path  App 

Ef  Total  App 

Bounces 

[W/m2] 

[W/m2] 

[W/m2] 

[W/m2] 

MWIR 

10 

7.13  x  10"3 

5.08  x  10"6 

2.19  x  HT6 

7.12  x  10“3 

LWIR 

10 

4.18  x  IQ"2 

5.90  x  10"4 

1.09  x  10"4 

4.13  x  10"2 

4-3.3  Position  2.  The  second  position  provides  a  side-on  view  of  the  air¬ 
craft,  this  time  at  a  closer  range,  509  m  compared  to  1, 119  m  for  the  first  position. 
Figures  4.19  and  4.20  show  the  view  of  the  aircraft  at  Position  2  from  the  observer’s 
position  in  the  MWIR  and  LWIR  bands,  respectively.  Again,  note  that  in  the  LWIR 
case  the  terrain  reflection  is  increased.  Also  note,  that  from  this  angle  (Position  2),  all 
four  exhaust  plumes  are  visible.  However,  as  the  hot  engine  internals  are  not  visible 
at  this  angle,  the  intensity  emitted  from  the  aircraft  is  much  less,  1852  W/sr  versus 
3800  W/sr  for  Position  1  in  the  MWIR  case.  However,  as  the  range  is  reduced,  the  ir¬ 
radiance  at  the  detector  increases  to  be  comparable  to  the  first  position,  a  coincidence 
of  the  two  ranges.  The  results  for  this  second  position  are  listed  in  Table  4.2. 

4.3.4  Position  3.  The  third  and  final  case  simulated  is  Position  3.  The 

aircraft  is  3  km  out  from  the  runway  at  a  height  of  150  m,  and  the  detector  is  viewing 
the  aircraft  from  front  on.  Figures  4.21  and  4.22  show  the  attitude  of  the  aircraft 
as  viewed  from  the  observer.  Note  the  large  contribution  from  the  background  and 
reduced  contribution  from  the  engines,  because  of  the  aircraft  orientation.  Table  4.3 
lists  the  results  from  the  simulations.  The  LWIR  case  was  simulated  several  times, 
whilst  varying  the  number  of  bounces,  similar  to  the  MWIR  simulation  at  Position 
1.  This  time,  the  integrated  irradiance  does  vary  with  the  number  of  bounces.  The 
observation  angle  in  this  case  can  support  multiple  bounces  and  the  observed  irra¬ 
diance  increases  as  each  bounce  is  included,  up  to  four  bounces  where  all  paths  of 
appreciable  radiance  exit  to  the  background.  These  multiple  bounce  paths  occur  as 
the  reflected  rays  bounce  between  the  fuselage,  wings  and  engines.  Figure  4.23  shows 
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(b)  Aircraft  plus  Background 

Figure  4.19:  The  two  MWIR  images  show  (a)  the  background  removed,  and  (b)  the 
background  included  to  show  contrast  between  the  aircraft  emission  and  reflection, 
and  the  background  irradiance  at  Position  2.  Colormap  units  are  [W/(m2  —  sr)]. 
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(a)  Aircraft  Only 


(b)  Aircraft  plus  Background 

Figure  4.20:  The  two  LWIR  images  show  (a)  the  background  removed,  and  (b)  the 
background  included  to  show  contrast  between  the  aircraft  emission  and  reflection, 
and  the  background  irradiance  at  Position  2.  Note  the  increased  contribution  in  the 
LWIR  band  from  the  ground  reflection  off  the  aircraft.  The  ground  (302. 5K)  has  its 
peak  emission  at  about  10/im,  in  the  middle  of  the  LWIR  band.  Colormap  units  are 
\W/(m2  —  sr)]. 
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Table  4.3:  Position  3  Detector  Irradiance. 


Band 

# 

Bounces 

Cf  Tgt  App 

[W/m2] 

Ee  Bg  App 

[W/m2] 

Cf  Path  App 

[W/m2] 

Cf  Total  App 

[W/m2] 

MWIR 

10 

5.73  x  10"4 

7.09  x  10"7 

4.03  x  10"7 

5.72  x  10"4 

LWIR 

1 

4.89  x  10"3 

8.25  x  10"5 

2.77  x  10"5 

4.83  x  10“3 

2 

5.27  x  10“3 

8.25  x  10“5 

2.77  x  10“5 

5.22  x  10“3 

3 

5.38  x  10"3 

8.25  x  10“5 

2.77  x  10"5 

5.33  x  10“3 

4 

5.43  x  10“3 

8.25  x  10“5 

2.77  x  10“5 

5.38  x  10“3 

10 

5.43  x  10“3 

8.25  x  10“5 

2.77  x  10“5 

5.38  x  10“3 

the  different  signature  between  a  single  bounce  and  ten  bounces,  where  the  single 
bounce  signature  is  10%  less  than  the  multiple  bounce  simulation. 

4-3.5  A340-300  Simulation  Summary.  The  simulations  which  calculate  the 
IR  signature  of  the  A340-300  aircraft  all  produce  similar  values  of  irradiance  at  the 
detector.  The  strongest  signal  is  from  Position  2  in  the  LWIR  band,  4.13  x  10-2  W/m2, 
and  the  smallest  is  from  Position  3  in  the  LWIR,  5.72  x  10-4  W/m2.  The  MWIR 
irradiance  was  on  the  order  of  10-3  W/m2.  Thus,  when  applying  the  SNR  condition, 
a  detector  must  have  an  NEI  no  greater  than  5  x  10-5  W/m2  to  detect  and  track  the 
aircraft  at  all  three  positions. 

4-4  Chapter  Summary 

This  Chapter  has  shown  the  results  of  the  cavity  blackbody  source,  IR  Range, 
and  the  A340-300  engagement  simulations.  The  cavity  blackbody  simulations  demon¬ 
strated  that  the  propagation  algorithm  is  numerically  accurate,  and  thus,  provides  a 
correct  result.  The  IR  Range  simulation  showed  that  the  propagation  algorithms, 
through  increasing  the  number  of  bounces,  were  able  to  more  accurately  predict  the 
IR  signature  of  the  cylinder  model  when  compared  to  the  LCAIR  algorithm  of  [3] 
and  [17].  Finally,  the  A340-300  simulations  show  the  utility  of  the  IR  signature  sim- 
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(a)  Aircraft  Only 


(b)  Aircraft  plus  Background 

Figure  4.21:  The  two  MWIR  images  show  (a)  the  background  removed,  and  (b)  the 
background  included  to  show  contrast  between  the  aircraft  emission  and  reflection, 
and  the  background  irradiance  at  Position  3.  Colormap  units  are  [W/(m2  —  sr)]. 
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(b)  Aircraft  plus  Background 

Figure  4.22:  The  two  LWIR  images  show  (a)  the  background  removed,  and  (b)  the 
background  included  to  show  contrast  between  the  aircraft  emission  and  reflection, 
and  the  background  irradiance  at  Position  3.  Note  the  increased  contribution  in  the 
LWIR  band  from  the  ground  reflection  off  the  aircraft.  The  ground  (302. 5K)  has  its 
peak  emission  at  about  10 fim.  Colormap  units  are  [I V/(m2  —  sr)]. 
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(a)  Single  Bounce 
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(b)  10  Bounces 

Figure  4.23:  The  two  LWIR  images  show  the  increased  signature  level  with  the 

multiple  bounce  simulation  at  Position  3.  Note  the  reduced  number  of  black  facets 
and  the  specular  reflection  off  the  fuselage  in  (b).  Colormap  units  are  [W/(m2  —  sr)]. 
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ulation  code,  by  calculating  the  irradiance  of  the  aircraft  at  the  optic  of  the  observer, 
for  several  scenarios. 

These  results  will  be  analyzed  and  discussed  in  detail  in  Chapter  V. 
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V.  Conclusions  and  Recommendations 


The  aim  of  this  thesis  was  to  show  how  BRDFs  affect  IR  signature  models.  Chap¬ 
ter  IV  showed  the  effects  of  including  the  BRDF  versus  a  specular  reflection 
assumption.  For  the  most  part,  the  specular  assumption  worked  quite  well  for  the 
bare  aluminum  and  gloss  black  surfaces  as  their  reflectance  distributions  were  very 
close  to  specular.  Secondly,  the  effect  of  the  number  of  bounces  was  shown  and  proved 
important  when  a  highly  reflective  surface  was  simulated.  Recall  that  each  time  the 
reflected  radiance  term  dominates,  through  a  high  reflectance,  p(0),  value,  another 
bounce  must  be  calculated  to  ensure  the  observed  radiance  properly  evaluated. 

Additionally,  computation  time  was  investigated  and  was  highly  dependent 
upon  the  number  of  paths  carrying  appreciable  radiance.  For  a  perfectly  specu¬ 
lar  surface,  there  is  only  one  path  per  facet,  and  for  a  diffuse  surface,  every  path 
(#f  acetsA(#bounces))  is  important.  This  was  shown  in  the  cavity  blackbody  analysis, 
where  the  simulation  was  stopped  after  four  bounces. 

Considerable  work  was  also  performed  to  optimize  the  sampling  rate  to  ensure 
the  BRDF  was  sampled  appropriately,  ensuring  a  correct  result,  whilst  keeping  the 
simulation  time  to  a  minimum.  It  is  possible  to  write  an  algorithm  that  is  capable 
of  solving  every  combination  of  BRDF,  geometry  and  required  accuracy.  However,  it 
would  not  be  efficient,  but  with  a  few  simplifications  (similar  to  the  ones  described  in 
Chapter  III),  the  algorithm  could  be  tailored  to  provide  a  much  improved  run  time. 

For  bare  aluminum  and  gloss  black  painted  surfaces,  the  specular  assumption 
works  well  and  negates  having  to  integrate  the  full  BRDF. 

5.1  LCAIR  Algorithm  Summary 

Overall,  LCAIR,  with  its  single  bounce  algorithm  and  purely  specular  assump¬ 
tion,  predicted  the  observed  radiance  well  for  situations  where  the  following  two  con¬ 
ditions  are  satisfied. 
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•  The  specular  bounce  from  the  first  facet  must  hit  a  facet  with  a  high  emissivity, 
whether  background  or  another  model  facet. 

•  The  reflectance  distribution  of  the  surface  must  satisfy  the  specular  assump¬ 
tion.  This  means  that  the  BRDF  of  the  surface  must  resemble  a  delta  function 
centered  on  the  specularly  reflected  direction. 

When  these  two  conditions  are  met,  then  the  LCAIR  algorithm  calculates  the  cor¬ 
rect  result  rather  well.  However,  when  one  or  both  of  the  conditions  are  not  met, 
then  the  result  is  incorrect.  The  second  requirement  was  not  demonstrated,  however, 
Sections  3.3  and  3.3.3  have  dealt  with  the  form  of  the  BRDF,  regarding  computation 
time,  sampling  resolution  and  hemispherical  irradiance  integral  approximations. 

A  simple  improvement  to  the  LCAIR  algorithm  would  be  to  include  multiple 
bounces.  The  specular  assumption  removes  the  BRDF  integration  overhead  by  simply 
using  reflectance  data  on  each  bounce.  Increasing  the  number  of  bounces  from  one 
to  many  would  not  affect  the  run  time  appreciably  as  most  facets  only  support  a  few 
bounces.  This  would  still  allow  trend  analysis  as  a  function  of  azimuth  (that  was 
demonstrated  in  [17])  to  be  performed  in  a  comparable  timeframe. 

5.2  Lessons  Learned 

Similar  to  the  lack  of  a  priori  information  that  hinders  an  algorithm’s  perfor¬ 
mance,  the  same  can  be  said  of  the  education  process,  and  hindsight  is  always  20/20. 

Most  open  source  or  even  commercial  ray-tracing  programs  are  written  in  C. 
The  decision  to  write  a  ray-tracer  in  Mat  lab®  proved  problematic,  and  more  time 
than  perhaps  was  warranted  was  spent  on  the  3D  wireframe  and  ray-tracing  portions 
of  the  software.  This  decision  was  made  due  to  the  author’s  unfamiliarity  with  C.  The 
extra  effort  involved  with  coding  in  Mat  lab®  compared  with  the  learning  curve  of  C 
meant  that  time  was  spent  on  the  framework  of  the  code,  rather  than  implementing 
the  propagation  algorithms  through  exercising  the  model.  An  improvement  to  this 
would  be  to  utilize  the  power  of  a  commercial  product  and  write  a  plug-in  (or  similar 
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code)  to  interface  with  the  3D  and  ray-tracing  engine.  Matlab®  is  not  optimal  for 
complex  3D  graphics  when  thousands  of  facets  are  involved.  A  move  to  C  or  another 
software  product  is  recommended. 

Additionally,  time  was  spent  attempting  use  the  radiosity  algorithm  because  of 
the  property  that  it  calculates  a  solution  for  every  observation  angle.  The  radiosity 
method  was  rejected  because  of  the  sampling  resolution  requirement  of  the  specular 
BRDFs.  A  better  approach  would  have  been  to  implement  the  POV  solution  first,  and 
then  experiment  later  with  global  illumination  algorithms.  Global  illumination  algo¬ 
rithms  typically  are  more  complex,  whether  numerical  quadrature  based  like  radiosity, 
or  probabilistically  based  like  the  Monte  Carlo  method. 

5.3  Recommendations  for  Future  Work 

The  first  recommendation  for  future  efforts  in  this  area  is  to  optimize  the  ray¬ 
tracing  algorithm,  which  will  provide  a  more  timely  solution.  OpenGL  software  rou¬ 
tines  provide  access  to  the  computer’s  graphics  hardware,  which  could  be  used  to 
ray-trace  the  scene  many  orders  of  time  faster.  One  method  previously  discussed  is 
to  write  a  plugin  (for  Blender3D  or  similar  product)  that  exploits  the  power  of  the 
3D  engine. 

Another  area  that  provides  promise  for  accurate  BRDF  calculations  is  a  tech¬ 
nique  called  Photon  Mapping.  It  is  a  probabilistic  technique  where  many  photons 
are  emitted  from  each  source  and  traced  through  the  model.  Each  surface  interaction 
(reflection)  is  recorded  in  a  photon  map,  and  the  photon  density  (observed  radiance) 
from  each  surface  is  estimated.  The  estimate  improves  as  the  number  of  photons 
increases.  This  method  can  provide  accurate  results,  however,  as  it  is  a  Monte  Carlo 
technique,  it  does  not  provide  the  quick  solution  required  for  trend  analysis. 

LCAIR  was  originally  developed  as  a  quick  IR  signature  trend  analysis  tool, 
to  be  used  instead  of  higher  fidelity  models  like  Spectral  and  Inband  Radiometric 
Imaging  of  Targets  and  Scenes  (SPIRITS).  However,  this  thesis  has  shown  that  having 
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to  integrate  the  BRDF  for  each  facet  is  time  consuming,  and  probably  unsuitable  for 
trend  analysis  work,  where  hundreds  of  observation  angles  are  required.  To  provide 
a  timely  trend  analysis  result,  the  specular  assumption  should  be  used,  however, 
this  assumption  requires  the  reflectance  distribution  of  the  material  to  be  highly 
specular.  This  requirement  limits  the  range  of  materials  that  can  be  simulated,  as 
not  all  materials  are  highly  specular. 

Aside  from  the  accuracy  of  the  algorithms,  the  LCAIR  software  has  not  been 
used  in  a  thorough  analysis  of  a  specific  aircraft’s  IR  signature.  Simulations  with 
physically  measured  surface  reflectance  distributions  and  correctly  measured  temper¬ 
atures  could  be  performed  to  analyze,  in  detail,  the  signature  from  an  aircraft  or 
other  targets  of  interest.  Other  additions  to  the  model  could  include  the  detector  sys¬ 
tem,  including  material  spectral  efficiency,  optics  transmissivity,  and  electronic  noise 
sources.  This  would  provide  a  whole  electro-optical  system  simulation,  which  would 
be  able  to  simulate  engagement  scenarios  and  evaluate  various  detection  strategies. 
Additionally,  including  an  imaging  detector  model,  whether  a  scanning  mirror  or  focal 
plane  array,  would  provide  even  greater  scope  for  system  simulation  and  analysis. 

5.4  Conclusion 

Several  algorithms  have  been  presented  and  implemented  for  simulating  IR  sig¬ 
natures.  Assumptions  of  specular  reflectance  can  reduce  computation  times  signifi¬ 
cantly,  however,  the  BRDF  must  very  closely  resemble  a  delta-function  at  the  specu¬ 
larly  reflected  angle. 

Trend  analysis  requires  a  fast  algorithm,  and  works  well  with  the  specular  as¬ 
sumption,  however,  a  balance  between  processing  time  and  accuracy  is  required. 

Finally,  multiple  bounce  calculations  are  important  when  the  geometry  is  en¬ 
closed  or  complex  (which  will  support  the  multiple  bounces),  and  also  when  highly 
reflective  surfaces  are  used  which  require  the  reflected  radiance  term  to  be  calculated 
with  additional  bounces. 
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Appendix  A.  Simple  Geometric  Analytical  Proofs 


Two  of  the  simulations  performed  to  provide  a  qualitative  check  of  the  algorithms 
are  included  in  this  appendix  for  reference.  They  were  chosen  as  an  analytical 
solution  can  be  gained  to  compare  with  the  numerical  result.  In  both  cases,  any 
contribution  from  the  background  is  ignored. 

A .  1  Two  Plates 

The  first  simulation  is  two  20  cm  x  10  cm  rectangular  plates  that  are  joined  at 
their  short  edges  at  right  angles.  Figure  A.l  is  a  depiction  of  the  setup.  The  first 
plate  has  a  temperature  of  800  K  and  is  a  perfect  blackbody  Lambertian  emitter. 
The  second  plate  is  a  perfect  diffuse  reflector  with  p  =  1.  The  temperature  is  300K, 
although,  it  does  not  matter  as  e  =  0.  The  observer  is  located  normal  to  the  reflecting 
surface  such  that  the  emitting  surface  is  not  visible. 

To  analytically  calculate  the  intensity  viewed  by  the  observer,  the  irradiance 
onto  the  second  plate  must  be  calculated.  Small  angle  approximations  cannot  be  used 
as  the  geometry  is  too  small  and  a  diffuse  reflectance  is  being  used,  thus  the  integral 
must  be  evaluated.  The  irradiance  on  the  second  plate  is  calculated  in  Equation  (A.l), 
and  the  radiance  emitted  from  a  perfectly  diffuse  surface  is  simply  the  irradiance 
multiplied  by  the  reflectance  and  divided  by  7r  (Equation  (A. 2)).  The  intensity  as 
viewed  at  the  observer’s  location  is  calculated  with  Equation  (A. 3)  from  the  emitted 
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T=300K 
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Figure  A.l:  Geometrical  setup  for  the  two  plates  simulation. 

The  observer  cannot  see  the  emitting  plate,  only  the  diffuse  re¬ 
flection  from  the  other. 

The  intensity  at  the  observer  is  then  21.53  W/sr.  The  results  from  the  simula¬ 
tion  is  shown  in  Figure  A. 2,  where  the  intensity  from  each  differential  piece  of  area  is 
shown.  The  output  from  the  simulation  is  22.05  W/sr  which  is  102%  of  the  analytical 
result. 


Figure  A. 2:  The  result  of  the  two  plates  simulation  is  shown. 
Note  the  gradient  profile  in  the  irradiance  on  the  second  plate. 
The  observer  is  looking  vertically  down  in  this  figure.  Colormap 
units  are  \W / sr\. 
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Figure  A. 3:  The  geometry  of  the  second  simulation  is  shown. 

The  observer  is  now  at  a  45°  angle  and  the  reflectance  is  now 
p  =  0.5. 

A. 2  Open  Box  with  Rear  Wall  Emitting 

The  second  analytical  simulation  starts  with  the  same  two  plates  as  before, 
however,  they  are  now  surrounded  by  a  box.  Figure  A. 3  shows  the  geometry  of  this 
situation.  The  emitting  plate  is  the  rear  wall  of  the  box  and  the  reflecting  plate  is  the 
blue  colored  wall.  The  observer  has  rotated  to  45°  from  normal.  The  side  walls  of  the 
box  shield  the  observer  from  the  direct  emission  of  the  emitting  wall,  thus,  the  only 
signal  that  the  observer  can  see  is  the  reflection  again,  but  this  time  the  reflectance 
of  the  second  plate  is  p  =  0.5.  Thus,  the  observed  intensity  contains  reflected  and 
self-emitted  radiance.  The  analytical  solution  is  shown  in  Equations  (A. 4)  and  (A. 5). 
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Figure  A. 4:  The  result  of  the  open  box  simulations  is  shown. 

In  this  simulation,  the  observed  intensity  is  a  combination  of 
reflected  and  self-emitted  radiance.  Colormap  units  are  \W / sr]. 

The  analytical  solution  is  13.7  W/sr  and  the  computer  simulation  using  the 
radiosity  algorithm  is  13.1  W/sr.  Again,  the  simulation  is  accurate,  96%  in  this  case. 
To  gain  greater  accuracy,  the  number  of  differential  areas  must  be  increased.  In  this 
simulation,  only  100  areas  are  calculated. 
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Appendix  B.  Blackbody  Cavity  Source  Calculations 


This  appendix  details  the  analytical  calculations  necessary  to  produce  Figure  3.33 
during  the  blackbody  cavity  source  simulation. 


B.l  s/S  Ratios 

Table  B.l  calculates  the  s/S  ratios  for  input  into  Equation  (3.12). 


B.2  Emissivity  Calculations 

An  example  calculation  of  the  analytical  emissivity  and  simulated  emissivity, 
for  a  specific  s/S  ratio  is  included  below  for  clarity.  To  start,  the  cavity  with  four 
facets  removed  is  chosen  for  the  simulation.  The  internal  emissivity  is  £  =  0.5,  and 
the  equations  are  taken  from  Section  4.1.  The  internal  temperature  of  the  cavity  is 
800  Kelvin. 


s 

S 

S0 

k 


£o 


0.0487[m2] 

3.106[m2] 


47rr2  =  47t(0.5)2  =  7r[m2] 
s  s 
S~S~o 


(1-s) 


(1-0.5) 


0.0487  0.0487 


3.106 
0.5 


7T 


=  8.88  x  10 


-5 


e[1_f]+(f)  0.5  [1-0.0157] +  0.0157 

£g(l  +  k)  =  0.9845(1  +  8.88  x  10"5)  =  0.9846 


=  0.9845 


Analytically,  the  emissivity  of  the  cavity  should  be  0.9846.  After  the  simulations  were 
performed,  the  emissivity  of  the  cavity  was  calculated  using  the  observed  radiance 
and  the  analytical  radiance  of  the  aperture.  The  cavity  was  set  to  a  temperature 
of  800K.  The  following  calculates  the  intensity  that  should  be  visible  from  a  perfect 
blackbody  with  an  aperture,  s  =  0.0487 m2 . 
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Figure  B.l:  Eight  different  aperture  sizes  were  constructed  for  the  simulation,  and 
are  displayed  above.  Each  figure  shows  the  size  of  the  aperture  in  relation  to  the 
radius  of  the  cavity,  the  s/S  ratio.  The  sub-caption  for  each  figure  lists  how  many 
facets  were  removed  to  create  the  aperture.  The  line  from  the  center  of  the  cavity  is 
the  view  direction  vector,  it  starts  at  the  center  of  the  cavity  and  extends  out  through 
the  center  of  the  aperture. 
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Table  B.l:  Gouffe  variables  for  each  of  the  Black- 

body  Cavity  Simulations. 
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A  perfect  blackbody  should  be  emitting  360.04  [IT/.sr]  on  the  centerline.  From  above, 
the  direct  emission  of  a  e  =  0.5  greybody  should  be  180.02  [IF/sr],  which  is  the 
simulation  result  after  zero  bounces.  After  one  bounce  the  viewed  intensity  from 
the  aperture  is  270.54  [IF/sr],  The  ratio  of  the  viewed  intensity  and  the  analytical 
blackbody  intensity  is  the  apparent  emissivity  of  the  cavity,  in  this  case  0.75.  This 
data  point  is  the  second  point  on  the  single  bounce  curve  of  Figure  3.33. 
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To  complete  the  rest  of  the  simulations  to  £11  in  Figure  3.33,  this  process  is 
repeated  for  all  combinations  of  the  eight  s/S  ratios  and  the  four  different  number  of 
bounces. 
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Appendix  C.  Ray  Tracing 


This  appendix  demonstrates  the  technique  of  ray  tracing ,  and  how  it  is  used  to 
sample  the  hemispherical  view  of  each  facet.  The  ray  tracing  technique  is  taken 
from  [9,15].  These  techniques  are  the  basis  of  classical  computer  graphics  algorithms, 
which  are  used  to  render  artificial  scenes  from  a  3D  wireframe. 

The  propagation  algorithm  needs  to  know  what  model  facets  are  viewable  from 
each  of  the  facets  in  the  wireframe  model.  To  achieve  this,  the  2D  hemisphere  is 
sampled  to  produce  a  hemispherical  view  matrix  for  each  facet.  The  matrix  has  a 
sample  every  degree  in  9  and  </>,  making  360  x  90  =  32400  samples  per  facet. 

C.l  Wireframe  Model  Primitives 

The  wireframe  model  contains  the  vertex  and  normal  vector  information  for 
every  facet.  Vertices  must  be  coplanar,  which  produces  a  flat  surface  which  creates 
either  a  triangle  or  quadrilateral.  The  ray  tracing  algorithm  must  then  be  able  to 
calculate  intersections  with  triangles  and  quadrilaterals. 

C.2  Ray  Tracing  Algorithm 

To  calculate  what  facet  is  visible  in  a  certain  direction  from  a  facet,  the  algorithm 
must  know  the  starting  position  and  direction  of  the  sample  ray.  Figure  C.l  defines 
the  geometry  that  will  be  assumed  for  the  description  of  the  ray  tracing  routine. 
Three  sample  rays  are  also  defined  in  the  figure  to  aid  in  the  description. 

Once  a  sample  direction  and  the  starting  location  (facet  A’s  center)  are  selected, 
every  other  facet  must  be  tested  to  see  if  that  vector  passes  through  the  facet.  To 
do  this,  the  algorithm  calculates  the  intersection  between  the  sample  vector  (ray) 
and  the  infinite  plane  that  the  four  vertices  create.  Then,  the  intersection  point  is 
checked  to  see  if  it  is  inside  the  polygon  or  not.  If  so,  the  distance  from  facet  A’s 
center  to  the  intersected  polygon’s  center  is  recorded.  The  next  polygon  is  checked 
against  the  sample  ray,  and  the  process  continues  until  all  polygons  are  checked  for 
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3 


facet  A 


Figure  C.l:  Three  sample  rays  are  drawn  to  aid  in  the  ray 

tracing  example.  The  facets  found  by  the  rays  are  facet  B,  back¬ 
ground  and  facet  C  for  sample  rays  1,  2  and  3,  respectively. 

an  intersection.  The  polygon  that  has  the  closest  intersection  to  the  parent  facet  is 
then  the  facet  that  is  in  view. 

For  example,  sample  ray  #1  will  find  facet  B.  Sample  ray  ^2  will  not  find  any 
of  the  model  facets,  and  thus,  background  will  be  assigned  in  that  direction.  Sample 
ray  #3  is  slightly  more  complex  as  the  piecewise  algorithm  will  find  facets  C  and 
D.  Clearly  facet  C  is  the  closest  and  should  be  assigned  to  this  ray.  The  algorithm 
calculates  the  distance  to  each  facet  and  chooses  the  closest  one,  facet  C. 

This  process  is  repeated  for  each  facet  in  the  model  and  creates  a  data  structure 
where  the  hemispherical  view  for  each  facet  as  a  function  of  6  and  <fi  is  indexed 
against  the  other  facets.  The  variable  names  are  geomFile  for  the  model  facets  and 
bgGeomFile  for  the  background  facets.  This  precalculation  creates  a  simple  lookup 
table  that  the  algorithms  evaluate  when  sampling  or  integrating  the  hemisphere  to 
allow  the  simulation  to  render  a  solution  faster.  This  process  is  the  most  demanding 
computational  process  in  the  simulation. 
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